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ANALYSIS  OF  SURFACE  ABLATION  OF  NONCHAERING  MATERIALS 

WITH  DESCRIPTION  OF  ASSOCIATED  COMPUTING  PROGRAM 

'  By  Fred  W.  Matting  and  Dean  R.  Chapman 
Ames  Research  Center 


SUMMARY 


[a  generalized  method  is  presented  for  solving  the  prohlem  of  stagnation- 
point  heat  transfer  and  material  response^^for  hlunt  Bodies  experiencing  melt¬ 
ing  and  vaporizing  or  suhliming  ahlationTj  The  analysis  is  applicable  to 
wind-tunnel  and  flight  conditions  (with  body  forces  taken  into  account) ;  inter¬ 
nal  radiation  can  be  considered  or  the  material  can  be  assumed  opaque;  the 
analysis  can  be  used  for  different  planets.  During  entry  flights,  a  body  will 
start  in  the  free -molecule  regime,  pass  through  a  transitional  regime,  and 
finish  in  the  continuum  regime  of  gas  dynamics.  Approximate  equations, 
rationally  obtained,  are  presented  which  provide  "bridges"  between  the  free- 
molecule  and  continuum  regimes  for  such  quantities  as  the  convective  heating 
rate,  surface  shear,  heat  blockage,  and  mass  loss.  Several  Illustrative  exam¬ 
ples  (including  surface  chemical  reaction  cases)  show  the  applications  of  the 
analysis.  Comparisons  with  experiment  are  made  where  possible.  The  analysis 
has  been  machine  programmed  for  numerical  solutions  using  a  finite  difference 
scheme,  and  a  running  energy  balance  is  kept  as  a  check  on  accuracy^  Instruc¬ 
tions  are  provided  for  using  the  computing  program. 


INTRODUCTION 


The  development  of  heat  shields  for  space  vehicles  and  long-range 
missiles  has  motivated  a  marked  increase  in  the  effort  to  understand  the  pro¬ 
cess  of  ablation.  This  intensified  study  is  contributing  to  the  understanding 
of  natural  ablative  phenomena  that  occur  when  extraterrestrial  bodies  enter 
the  Earth’s  atmosphere.  Ablation  data  obtained  in  the  laboratory,  for  example, 
in  arc -jet  wind  tunnels,  do  not  duplicate  in  any  single  experiment  all  the 
conditions  of  entry  flights;  hence,  there  is  need  for  analytical  methods  of 
predicting  and  explaining  the  ablative  phenomena.  Before  such  methods  can  be 
applied,  however,  their  validity  and  accuracy  must  be  determined  by  comparing 
calculated  results  with  results  from  wind-tunnel  tests,  with  flight  data,  or 
with  post-flight  observations  of  a  man-made  or  natural  object  when  one  can  be 
recovered. 

A  generalized  method  is  presented  here  for  solving  the  combined  problem 
of  heat  transfer  and  material  response  for  the  stagnation  region  of  blunt 
bodies  experiencing  melting  and  vaporizing  or  subliming  ablation.  An  attempt 
has  been  made  to  describe  the  problem  mathematically  as  completely  as  possible 
in  order  to  obtain  nearly  exact  solutions.  This  required  that  the  analysis  be 
machine  programmed  for  numerical  solutions.  This  program  in  various  stages  of 


development  has  been  used  successfully  at  Ames  Research  Center  during  the 
past  three  years.  It  has  previously  been  employed  in  the  analysis  of  tektite 
ablation  (refs.  1,  2). 

A  number  of  the  equations  used  in  the  analysis  are  already  well 
established  (refs.  3-6) j  however,  some  of  the  equations  and  features  of  the 
analysis  are  new.  The  analysis  takes  account  of  the  fact  that  entry  bodies 
initially  fly  in  the  free -molecule  regime,  then  in  a  transitional  regime,  and 
finally  in  the  continuum  regime  of  gas  dynamics.  The  analysis  contains  for¬ 
mulas  for  the  transitional  regime  to  bridge  between  the  free -molecule  and 
contlnuim  regimes;  these  formulas  have  been  rationally  derived  from  simple 
models  and  are  believed  to  fill  an  important  gap  in  previous  analyses  of  small 
objects  entering  a  planetary  atmosphere. 

Several  options  in  the  analysis  and  associated  computing  program  are 
available.  Internal  radiation  in  the  body  is  accounted  for,  or  the  body  can 
be  assumed  to  be  opaque.  Flight  cases  as  well  as  wind-tunnel  cases  can  be 
calculated;  the  flight  cases  can  be  applied  to  any  planet,  provided  certain 
characteristics  of  the  atmosphere  are  known.  The  rear  boundary  conditions  for 
the  ablating  material  can  be  those  for  a  heat  shield,  or  the  aerodynamic  base 
heating  for  an  object,  such  as  a  tektite,  can  be  accounted  for.  The  ablating 
material  can  be  a  type  that  melts  and/or  vaporizes,  sublimes,  or  undergoes  a 
surface  chemical  reaction  in  the  ablation  process.  The  various  material  prop¬ 
erties  and  the  external  flow  conditions  can  be  put  into  the  computing  program 
arbitrarily  so  that  a  variety  of  ablation  research  problems  can  be  studied. 

To  illustrate  the  types  of  problems  that  can  be  handled  by  the  analysis 
(and  to  elucidate  some  of  the  significant  ablative  phenomena) ,  several  typical 
examples  are  presented.  These  are  calculated  by  the  numerical  computing  pro¬ 
gram  associated  with  the  analysis.  The  principal  features  of  the  numerical 
computing  program  are  described,  and  instructions  are  given  for  its  use. 


ANALYSIS  AND  METHOD  OF  SOLUTION 


In  this  section,  the  principal  emphasis  is  on  the  method  of  solution 
used  in  the  associated  computer  program.  Equations  are  presented  in  specific 
form  with  numerical  values  inserted  for  universal  constants.  In  performing 
calculations  for  various  materials  and  situations,  the  user  has  a  number  of 
parameters  (constants)  at  his  disposal;  typical  values  of  these  are  listed  in 
the  description  as  they  appear. 


Basic  Approach  and  Approximations 

The  analysis  is  concerned  with  the  problem  of  surface  ablation  near  the 
stagnation  point  for  transparent  or  opaque  materials  of  finite  thickness.  The 
analysis  is  thus  restricted  to  materials  that  undergo  surface  ablation  includ¬ 
ing  melting,  evaporation,  sublimation,  or  surface  chemical  reactions,  such  as 
depoljrmerization.  Chemical  reactions  involving  the  external  gases  are  not 
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considered.  The  houndary-layer  equations^  as  such,  are  not  solved,  hut 
results  of  known  solutions  are  used  to  obtain  heating  rates,  surface  shear, 
and  the  effects  of  mass  transfer.  Front  and  hack  surface  heating  rates  rmist 
he  solved  for,  in  order  to  determine  the  front  and  hack  face  boundary  condi¬ 
tions,  Two  types  of  heating  environments  are  treated.  One  is  of  constant 
velocity  and  stream  density,  such  as  generally  exists  in  an  arc -heated  wind 
tunnel;  the  other  is  of  time-varying  velocity  and  density  over  an  entry  flight 
(using  a  constant  lift/drag  ratio  hut  a  variable  mass).  In  this  latter  envi¬ 
ronment,  the  equations  of  motion  are  solved  simultaneously  with  the  heat- 
transfer  and  ablation  equations.  In  the  equations  of  motion,  the  ratio  of 
body  mass  to  the  product  of  drag  coefficient  by  frontal  area  appears,  and  this 
is  related  (empirically)  to  the  surface  recession. 


Conservation  Equations 


The  analysis  is  essentially  a  time -dependent  energy  balance  along  the 
stagnation  center  line  of  an  axisymmetric  blunt  body.  The  basic  equations  to 
be  solved  are  the  conservation  equations  for  energy,  mass,  and  momentum  for 
the  material  of  the  body,  written  in  a  simplified  form  that  is  valid  in  the 
stagnation  region.  The  energy  and  momentum  equations  have  been  further  sim¬ 
plified  by  neglecting  inertial  terms;  this  procedure  is  a  valid  approximation 
for  viscous  ablative  materials,  such  as  glass,  stone,  or  any  subliming  mate¬ 
rial.  The  curvilinear  coordinate  system  used  is  shown  in  sketch  (a).  The 

conservation  equations  are: 

Energy 

At  ^  -  St\  a  A  St  .  X 

vs  *  L  57  ■  J 

(1) 

where  F(y,t)  is  the  internal  radiation 
flux  term. 

Sketch  (a) 

Continuity  (for  constant  density) 
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Most  previous  analyses  have  taken  account  of  internal  radiation  or  vehi¬ 
cle  acceleration,  but  not  both.  The  acceleration  term  in  the  momentum  equa¬ 
tion  (3)  can  be  important  in  determining  the  rate  of  removal  of  melted 
material;  this  term,  as  written,  is  valid  in  the  stagnation  region. 


The  continuity  and  momentum  equations  can  be  put  into  quadrature  form. 

We  make  use  of  the  fact  that  _u  varies  linearly  with  x  near  the  stagnation 
center  line^  and  we  define:  u  =  (Su/Sx)jj_q.  Then^  near  the  center  line 


u  =  ux 


The  continuity  equation  in  quadrature  form  is: 


V  =  Vv,  -  2  /  u  dy^  (5 

^o 

The  value  of  v^  depends  on  surface  temperature  and  external  conditions;  its 
evaluation  is  described  below  in  the  section,  Front  Face  Velocity.  The  sur¬ 
face  recession  rate  and  surface  recession  are,  respectively. 


vsr  =  VT  =  IvbfI  =  W(y-RTr^'^) 


VgFl^ti 


Pressure  Gradient 

The  solution  to  equation  (5)  can  be  £ubstituted  into  equation  (l),  but  to 
solve  equation  (5)  we  must  first  obtain  u(y,t)  from  the  momentum  equation. 

We  write,  for  the  pressure  near  the  center  line, 

P  =  p(0)  +  I  p"(0)x^  (7) 

and  equation  (3)  becomes,  after  substituting  equation  (4), 


|.(,|)  =p..(0)  -  f  =  -P'. 


where  p"  is  then  the  negative  of  the  second  derivative  of  pressure  with  a 
corrGcision  for  fliG  body  force  acfing  on  fhe  material.  We  can  infegrafe  egua- 
tion  (8)  twice  and  obtain: 


u(y,t)  =  P"  J' 


^BF  yidy- 


'^BF  dy. 


where  (dynes/cm^)  is  the  x  derivative  of  the  absolute  value  of  the 

surface  shear ^  (dT-^/dx)^^— q .  In  P”  the  quantity^  is  the  acceleration^ of  a 
body  in  flight  (positive  for  increasing  speed) .  For  wind-tunnel  calculations 


where  K  is  the  fractional  time  that  the  center  point  is  initially  exposed  to 
approximately  stagnation  conditions,  and  and  are  values  selected  to 

give  a  realistic  damping  to  the  oscillation  during  entry.  Typical  values  for 
the  constants  that  have  been  used  for  calculating  a  tum'bllng_tektite  entry  are: 
iC  =  0.25;  =  R/10  cm;  Eis  =  1.0.  For  a  nontumbling  body,  K  is  unity  and 

the  other  constants  have  no  influence. 


Surface  Convective  Heat  Transfer 

In  order  to  evaluate  the  x  gradient  of  wall  shear,  ,  we  require 
first  an  evaluation  of  the  surface  convective  heat  transfer;  we  also  need  this 
to  determine  the  surface  boundary  conditions.  Instead  of  enthalpy  (cal/g), 
we  use  an  "enthalpy  velocity"  (km/sec),  which  is  defined  as 

=  0.00836  hs  (13a) 

V  =  0.0915  «/h7  (1313) 
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For  the  external  gas  we  use  an  average  specific  heat, 

(14) 


In  laminar  continuum  flow  we  can  evaluate  the  surface  convective  heat  transfer 
as  (ref.  7)  (with  a  vorticity  correction) 

-  0.00836  CpT^.)  (1  +  (15) 

where  A4  is  a  constant  and  Cq  is  a  vorticity  correction  (generally  small) . 
The  results  of  equation  (15)  agree  well  with  existing  experimental  data  over 
an  extended  range  of  enthalpy  potentials.  The  value  of  A4  for  Earth  entries 
is  approximately  1.1;  in  wind-tunnel  tests,  A4  is  evaluated  with  a  calorim¬ 
eter.  With  a  blowing  correction  we  have 

lllfC  = 

where  is  evaluated  in  equation  (28) . 

For  high  altitude  flight  or  for  rarefied  wind-tunnel  conditions,  it  is 
necessary  to  consider  the  free-molecule  regime.  For  surface  convection  in 
the  free-molecule  regime  we  use  a  Newtonian  type  of  approximation  (ref.  8, 
pp.  395-403). 


(16) 


^FM  “  0.0836  CpT^.)  (17) 

In  the  transition  regime  between  free-molecule  and  continuum  flow,  the  con¬ 
vective  heat  transfer  will  have  a  value  bridged  between  the  evaluations  in 
equations  (16)  and  (17) •  This  has  been  derived  from  a  simple  kinetic  theory 
model  (see  appendix  C )  yielding  the  result 

V  '  (1  -  (IB) 

When  the  tumbling  correction  is  applied^  the  convective  heat  transfer  at  the 
front  wall  is 


*l\|fW^U 

Equations  (15)  to  (19)  are  used  to  determine  the  convective  heat  transfer 
to  the  front  face,  hut  it  is  also  of  interest  to  calculate  some  related  quan¬ 
tities.  When  no  material  is  being  lost  to  the  vapor  state,  xjr  =  1,  and 
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equation  (l8)  specializes  to 


^o  = 


1  -  e 


■'Ifm/  ‘Iqi 


and  equation  (19)  specializes  to 


<30  =  'loo^u 


When  material  is  lost  to  the  vapor  state ^  \|r  ^  1;  we  can  consider  that  ^ 
operates  on  qoc  not  on  qpj^.  Then  a  modified  ^  can  he  obtained  to 

operate  on  q^^.  We  define  the  modified  blowing  factor^  as: 


(22a) 


(22b) 


On  multiplying  both  sides  of  equation  (22b)  by  we  have  also 


‘lw=  'l^'lo 


When_we  substitute  equations  (l6),  (l8),  and  (20)  into  (22a),  the  evaluation 
of  ^  becomes 

„Y-.  -'ifmA'IoA 


A|f  1  -  e 


1  -  e 


-ifmY 


Equations  (20 )  to  (2^)  are  alternate  forms  entirely  equivalent  to  equa¬ 
tions  (l8)  and  (l9)  •  The  quantities  qoo^  although  not  needed  in 

determining  the  convective  heat  transfer^  are  computed  as  quantities  of 
interest  conceptually. 

The  bridging  relation  given  in  equation  (l8)  or  the  alternate  forms  in 
equations  (20)  and  (24)  will  automatically  take  account  of  changes  of  heating 
rates  as  a  body  flies  from  one  regime  into  another  and  will  place  control  in 
the  appropriate  regime.  Comparisons  with  available  measured  data  are  shown 
in  appendix  C. 

For  normal  (nonrarefied)  wind-tunnel  conditions^  the  bridging  relations 
given  above  are  not  needed;  also  there  is  no  tumbling,  so  we  have  simply: 
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Ktu  =  1 


(25a) 


(2513) 

■30  >  V  <==50 

4,  =  tOo  (25a) 

In  the  heat-transfer  relations  given  above  from  equation  (15)  to  equa¬ 
tion  (25)^  the  bloving  parameter^  is  needed.  To  evaluate  \lr ^  ve  need  to 
know  the  relationship  betveen  the  equilibrium  vapor  pressure  and  temperature 
for  the  ablating  material^  p^^  =  p^g(T^).  We  will  write 


where  p^  is  the  modified  equilibrium  vapor  pressure  when  the  equilibrium  is 
shifted  by  the  presence  of  other  gaseous  materials  in  the  boundary  layer.  An 
example  of  this  is  the  suppression  of  vaporization  of  silica  due  to  the  pres¬ 
ence  of  oxygen  in  the  atmosphere.  This  suppression  effect  for  silica  is 
analyzed  in  reference  9  which  an  analytic^  but  implicit^  expression  is 
obtained  for  the  modified  equilibrium  vapor  pressure.  The  use  of  the  exponent, 
E7,  is  an  empirical  accounting  for  this  effect  which  yields  values  within 
several  percent  of  the  values  obtained  by  the  rigorous  analysis  of  reference 
For  silica,  the  value  E7  =  1.4  has  been  used.  For  a  number  of  other  mate¬ 
rials,  E7  =  1,  or  the  modified  equilibrium  vapor  pressure  is  the  usual  equi¬ 
librium  value.  For  Pt2  atmospheres)  we  use  a  hypersonic  approximation 
(twice  the  free -stream  dynamic  pressure  with  a  correction  factor) : 


^  AiDV^ 

101.3 


(27) 


where  a  typical  value  for  Ai  is  O.95,  and  the  number,  101. 3,  accounts  for 
the  units  in  the  equation.  Equation  (26)  is  a  limiting  form  for  continuum 
conditions;  P  could  be  based  on  the  actual  existing  vapor  pressure,  p.^, 
which  will  not  be  an  equilibrium  value,  but  in  the  limit  p^  approaches  p-y^ri 
from  below.  However,  p^  cannot  exceed  p-^^,  so  for  this  equation,  whatever 
value  p^  may  actually  take,  it  is  not  allowed  to  exceed  pt^ 

tion  (26;  .  This  has  been  done  by  giving  P  a  lower  limit  which  is  a  small 

positive  number  (lO”®).  So  P  is  defined  by  equation  (26)  down  to  its  lower 

limit.  The  quantities,  p^,  Py^,  Py^u?  considered  to  be 

evaluated  at  the  liquid  or  solid  surface. 


We  .can  represent  for  an  evaporation  or  sublimation  process  as 


+  E35 


(28a) 


P 


This  relation,  with  the  asymptote,  E35  =  O.O6,  gives  a  good  fit  to  a  number  of 
boundary -layer  solutions  (refs.  10-12).  Equation  (28a)  is  also  a  limiting 
form  for  continuum  conditions  since  P,  as  evaluated  in  equation  (26),  is  not 
based  on  the  actual  vapor  pressure,  p^.  The  quantity  l/P  is  a  function  of 
the  mass  loss  rate  due  to  vaporization  (see  eqs.  (38),  (39a),  and  (.40)  below). 
For  a  surface  chemical  reaction  we  use  the  form 


^  ^  -  +  E35  (28h) 

Bii 

P[1  +  (v^cAwFm)J 

where  v™  and  Vypjq  are  shown  evaluated  velow  in  equations  (38)  and  (39e). 

The  correction  using  v^c  and  takes  account  of  the  possible  reaction 

rate  control  of  the  mass  loss  rate  which  does  not  depend  on  the  (modified) 
equilibrium  vapor  pressure  (eqs.  (39^)  and  (40)).  The  quantity^  Bn  in  equa¬ 
tion  (28)  (0.95-1.55  used),  depends  on  the  ratio  of  molecular  weights  of  exter¬ 
nal  gas  to  blowing  gas  and  should  preferably  be  determined  by  experiment.  In 
the  absence  of  experiment,  we  can  estimate  Bn  as  follows.  We  can  define 


and  we  have 


constant 

Bn  - 

M 

where  the  constant  ^  0.7  to  0.8  and  n  2/3  to  3/4* 


(29) 


(30) 


Wall  Shear  Gradient 

We  are  now  in  position  to  evaluate  the  x  gradient  of  the  wall  shear, 

T  »  =  (dTv/dx)x=o-  consider  the  case  of  continuum  flow  with  no  blow- 

Ing  to  evaluate  Tq'  =  (dTo/dx)x=o-  Using  a  modified  form  of  Reynolds 
analogy,  we  can  write: 


I0C  Eig-To* 


(31) 


where  Kg  is  a  constant  that  depends  on  the  Prandtl  number  and  is  unity  when 
the  Prandtl  number  is  unity.  Using  equation  (ll)  we  obtain 


n/  2A2  _  Vq  ^Rn/^pi  (4^) 
==  “  lo^qoc^oo 


(32a) 
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As^cVoqIO- 

R(Ah)'s/^ 


(3213) 


where  A3  ^  1*45  is  a  typical  value.  We  apply  the  tumbling  or  oscillation 
correction,  (^1*  (l2))  to  and  obtain 


Tqc'  =  ""o'Ktu 


For  the  effect  of  blowing  on  wall  shear  in  continuum  flow  we  use  (ref.  11) 


The  value  of  Es  should  be  small;  unless  experimental  evidence  is  available, 
it  is  suggested  that  the  value  zero  be  used.  (Reasonable  answers  have  also 
been  obtained  with  Es  =  O.3  Bn.)  In  the  numerical  computing  program,  the 
value  of  P  that  is  inserted  into  equation  (3^)  (only)  is  arbitrarily  pre¬ 
vented  from  becoming  less  than  Ea  so  that  the  shear  blowing  factor  cannot 
become  unrealistically  large  when  the  actual  P  is  very  small.  Using  equa¬ 
tions  (32b,  33^  and  3^)^  w-e  obtain  our  expression  for  in  continuiom 

flow. 

,  _  836A..^qop.VooKt.ii»[l  +  (Eh/p)^] 

T-^C  “  I - 

Rn/P21  (V^  -  0.00836  CpT^) 

For  the  x  gradient  of  shear  in  free  molecule  flow  (which  is  unaffected 
by  blowing),  we  have  (ref.  8,  pp.  395-^03) 

AcmWoo^lO^iCtu 


where  is  the  x  momentum  accommodation  coefficient.  For  the  bridging 

between  free-molecule  and  continuum  wall  shear,  we  use  essentially  the  same 
model  as  that  used  for  heat-transfer  bridging  (see  appendix  C). 


1  -  e 


“'^vFm/^wc^ 


The  evaluations  of  p"  and  given  above  along  with  knowledge  of  the 

temperature  distribution  enable  us  to  solve  equation  (9)  for  ff(y^t)  which 
can  be  substituted  into  equation  (5)- 
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Front  Face  Velocity 


To  solve  equation  (5)  we  need  to  evaluate  the  front  face  velocity^  Vy 
For  the  front  face  velocity  under  continuum  conditions^  we  use  the  so-called 
Lewis  analogy  (Le  =  l) ,  which  states  that  the  ratio  of  mass  diffusion  to  con 
centration  "gradient"  is  equal  to  the  ratio  of  continuum  heat  transfer  to 
enthalpy  potential: 


0.00836 

■pMP(V^  -  0.00836  CpT.^^) 


(38) 


An  equivalent  form  of  equation  (38)  is  given  as  equations  (31)  and  (32)  of 
reference  9  (see  also  ref.  I3) •  Equation  (38)  is  a  limitingform  for  contin¬ 
uum  conditions^  "because  of  the  way  P  is  evaluated  in  equation  (2  ).  (As 
noted  previously,  P  is  evaluated  at  the  solid  or  liquid  surface.)  Equa¬ 
tion  (38)  also  contains  empiricism  in  the  evaluation  of  t  (eq.  (28)).  In 
reality,  it  is  expected  that  the  diffusion  rate  should  reach  a  maximum  value 
with  a  very  small  P,  if  P  were  "based  on  the  actual  pressure  at  the  surface, 
p  .  This  would  also  require  that  the  V  asymptote  "be  zero.  With  a  finite 
asymptote  (which  seems  to  fit  existing  data),  the  calculated  diffusion^ 
rate  becomes  unrealistically  large  for  small  P.  Under  these  conditions  it 
can  he  expected  that  free -molecule  or  reaction-rate  control  will  generally 
prevail  (see  eq.  (40)),  so  that  an  inaccurate  calculation  of  the  diffusion 
rate  for  these  conditions  will  have  little  effect  on  the  net  rate  calculated. 


For  the  calculation  of  the  front  face  velocity  in  the  free-molecule  or 
rate -controlled  regime,  we  distinguish  two  cases:  (l)  evaporation  or  sublima¬ 
tion,  and  (2)  a  chemical  reaction  such  as  a  depolymerization.  For  the  evapo¬ 
ration  or  sublimation  case,  we  have,  from  kinetic  theory,  the  Langmuir 
equation  (Cl) .  With  constants  evaluated  to  account  for  our  units  we  have: 


-  I  _  ^^•3Ap.vPve 

'%FMI  "  p  ^  T^ 


(39a) 


Using  equation  (29)  we  write 

1%.fm1  =  239  (acv/^) 

We  now  use 


and  get 


(39b) 

(39c) 

(39d) 
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Equation  (39^)  is  the  form  we  use  in  our  calculations.  When  the  molecular 
weight  of  the  external  gas  is  29-1^  the  accommodation  coefficient  can  he 

used  directly  in  equation  (39^) ;  otherwise  the  accommodation  coefficient  must 
he  corrected  according  to  equation  (39c) . 

For  the  chemical  reaction  case_,  we  use  an  Arrhenius  rate  form  for  a  first 
order  reaction 


^wFmI  " 


“^tAw 


(39e) 


In  the  coordinate  system  used^  hoth  v^c  will  he  negative  quantities. 

The  bridging  equation  between  the  free -molecule  or  reaction-rate  controlled 
regime  and  the  continuum  or  diffusion  controlled  regime  turns  out  to  he  (see 
appendix  C): 


1 

VwFM 


+ 


1 

■^WC 


( toa) 


^wFM'^wc 


(40h) 


The  use  of  equation  (40)  automatically  places  the  front  face  velocity  in  the 
appropriate  controlling  regime:  the  diffusion-controlled^  the  rate -controlled^ 
or  the  transitional  regime. 


Internal  Radiation 

In  the  energy  equation  (l)^  we  require  the  evaluation  of  the  internal 
radiative  flux^  F(y^t).  In  the  evaluation  of  F^  we  assume  either  an  opaque 
body  (F  =  O)  or  a  transparent  gray  body  (two  shades  of  gray;  one  for  the 
incoming  gas  cap  radiation^  another  for  the  absorption-emission  of  the  heat- 
shield  material).  A  third  alternative  is  to  approximate  the  transparent  gray 
body  by  a  semitransparent  body  that  is  opaque  internally^  but  has  a  variable 
surface  emissivity. 

The  evaluation  of  the  radiative  flux  for  the  gray  transparent  body  is 
given  below;  it  is  similar  to  that  of  reference  l4  which  treats  scattered 
radiation  and  uses  an  exponential  attenuation.  The  present  evaluation  con¬ 
siders  one  reflection  from  the  front  and  the  rear  surfaces^  which  is  a  good 
approximation  for  materials  that  absorb  well. 
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F(y,t) 

2n^acT 


py  P^BF 

=  J  T^(Ti)E2[a(y  -  •n)]dTi  -  J  T^(Ti)E2[a(Ti  -  y)]dTi 


P^BF 

Reff  J  T4(Tl){E2[a(y  +  ri)]  -  E2[a(2ygp  -  (y  +  n))]}  <^1 


2n2aCT  ®  "  ^eff® 


^eff  effective  coefficient  of  reflection  for  planar  radiation  and  is 

related  to  the  maximum  emisslvlty  hy  the  relationship  (ref.  15) 

1  -  Heff  =  ^  (^: 

The  second-degree  exponential  integral  £2(2)^  used  as  the  attenuation  func¬ 
tion  in  the  integrals  of  equation  T^l)^  is  (ref.  l6^  appendix  l): 


/oo  g-zxi 

dxi  = 


In  actual  numerical  calculations,  equation  (4l)  is  used  directly  to  evaluate 
F  at  the  front  and  hack  faces  only.  The  quantity,  g  =  bF/by,  is  obtained  hy 
taking  the  derivative  of  equation  (4l),  and  g  is  numerically  evaluated  for 
use  in  equation  (l) .  Internal  values  of  F  are  then  obtained  hy  numerical 
quadrature  of  the  flux  derivative,  g.  The  gas  cap  radiative  flux,  q^,  is 
evaluated  with  an  empirical  approximation: 

qp  =  E4m)®^V®®iCt„  (44 


The  form  of  equation  (44)  is  deduced  from  experimental  correlations  presented 
in  reference  I7.  Input  constants  E4,  E5,  and  Eg  for  a  given  environment 
should  he  selected  to  fit  available  data.  The  level  of  radiation  and  the 
surface  reflectivity  are  both  accounted  for  in  the  evaluation  of  E4j  Eg  can 
vary  from  O.5  for  no nequl librium  radiation  to  I.7  for  equilibria  radiation 
in  air,  while  Eg  can  vary  from  5  fo  8  for  air.  After  the  exponents  are 
selected,  E4  should  be  chosen  to  give  the  proper  level  of  radiation.  Values 
of  the  constants  that  have  been  used  for  Earth  tektlte  entries  are: 


E4  =  0.76x10' 


0*5:  Eg  —  7 • ^ 


In  the  calculation  for  a  semitransparent  body,  the  material  is  considered 
to  be  internally  opaque,  and  the  front  surface  emisslvlty  is  varied  in  an 
appropriate  manner  with  the  thermal  thickness  of  the  temperature  profile. 

This  variation  is  derived  by  assuming  an  exponential  temperature  distribution 
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near  the  wall.  This  permits  a  closed  form  integration,  and  the  result  is 
further  approximated  to  the  following  simple  form; 


^FP  ” 


'max 


^max 


1  + 


2.k 

oA 


1  + 


E 


33 

A 


(i^5) 


where  E33  =  2.4/a  and  the  thermal  thickness,  A,  is  evaluated  in  equa¬ 
tion  (63).  (if  ibT/bj)y.  >  0  or  A  <  0,  then  epp  is  assigned  the  value  of 
^max’)  With  this  representation,  the  back-face  emissivlty  epp  is  assumed  to 
be  constant  and  F  =  0.  The  semitransparent  approximation  gives  virtually 
the  same  results  for  most  ablation  characteristics  as  the  transparent  case 
(see  table  l),  although  the  Internal  temperature  profiles  do  not  agree  closely. 
It  greatly  reduces  computing  machine  time,  however;  hence  it  is  used  for  most 
calculations  when  an  accuracy  of  the  order  of  10  percent  is  adequate. 


Boundary  Conditions 

The  evaluations  of  all  the  terms  in  the  energy  equation  (l)  have  been 
shown,  so  this  equation  is  in  a  form  to  be  solved  for  T(y,t) .  Boundary  con¬ 
ditions  are  needed  for  equation  (l),  and  these  are  determined  in  the  standard 
manner  by  writing  surface  energy  balances  for  the  front  and  back  surfaces, 
providing  for  the  appropriate  differences  between  the  opaque  and  transparent 
cases.  Options  are  provided  in  determining  the  rear  surface  boundary 
conditions  as  shown  below. 

The  energy  balance  for  the  front  surface  is  written: 


phvV^  +  q^j.  +  (1  -  Bi6)(qp 


(46) 


where 

Bis  =  0  for  the  opaque  and  semitransparent  cases 
Bi6  =  1  for  the  transparent  case 

In  the  coordinate  system  used^  the  front  face  remains  at  the  origin.  As 
material  is  lost^  the  location  of  the  hack  face  recedes  to  smaller  values  of 
j-QY’  region  in  which  we  are  solving  equation  (l)  becomes  smaller.  So^ 

to  know  the  location  of  the  hack  face^  at  any  time^  t^  the  system  of  equations 
must  have  heen  solved  up  to  time,  t. 

In  writing  the  hack  surface  energy  balance,  we  distinguish  between  a 
hack  surface  exposed  to  flight  conditions  and  a  hack  surface  in  contact  with 
a  hacking  material.  For  the  hack  surface  exposed  to  flight  conditions  (as, 
for  example,  with  a  tektite),  we  write  the  following  hack  surface  energy 
balance: 
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where  Bis  is  zero  or  unity  as  noted  above.  The  hack  surface  convective  heat 
transfer,  qgp,  maybe  a  function  of  surface  temperature.  We  relate  to 

the  front  surface  convective  heat  transfer  (without  blowing) ,  empirically 
as  follows: 


;bf _ 

-  I^Bf) 


%/f'^ 


1  +  (Ei9  “  i)(  1 


>DV  Ri 


where  Rb/F  is  defined  as  the  ratio  of  base  to  front -face  laminar  convective 
heating  when  normalized  by  the  respective  enthalpy  potentials.  To  account  for 
transition  to  turbulent  flow  at  the  base^  we  use: 


^EF  turbulent 
^BF  laminar 


(^9) 


Typical  values  of  constants  that  have  been  used  for  tektite  entries  are: 

%/F  0.01;  Ei3  =  0.01;  Eig  =5-0.  When  the  quantities,  Rg/P  and  6^^,,  are 

assigned  the  value  zero,  the  back  boundary  condition  is  adiabatic. 

For  numerical  computations  of  the  transparent  case,  the  computing  program 
has  been  arranged  so  that  it  is  possible  to  modify  equations  (1*6)  and  (I+7)  by 
assigning  some  radiant  energy  to  the  surface  energy  balances.  The  reasons  for 
this  are  explained  in  the  section,  METHOD  OF  THE  HUMEIRICAL  PROGRAM,  Boundary 
Conditions  for  the  Transparent  Case,  equations  (97)  (98) > 

The  other  back  bovindary  condition  of  interest  is  used  for  heat-shield 
calculations.  We  assxame  a  backing  material  which  acts  as  a  heat  sink  and  is 
at  uniform  (but  increasing)  temperature,  For  an  opaque  heat  shield  the 

heat  transfer  is  entirely  by  conduction  to  the  backing  material;  for  a  trans¬ 
parent  heat  shield,  we  assume  that,  in  addition  to  conduction,  the  radiative 
flux,  F(ygp.,t)  is  entirely  absorbed  by  the  backing  material.  Then,  in  place 
of  equations  (47)  to  (50)  we  write 


Bi6F(y3p,t) 


(51) 


where  c  is  the  heat  capacity  per  unit  area  of  the  backing  material.  In  the 
limiting  case  with  H  =  0,  the  back  boundary  condition  is  adiabatic,  and 
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equation  (51)  should  be  replaced  by  equations  (4?)  and  (48)  with  Rp/p  and 
€pp  assigned  the  value  zero. 

It  should  be  noted  that  the  surface  temperatures  generally  cannot  be 
specified  a  priori,  as  they  depend  partly  on  external  conditions  that  may 
change  with  time.  The  computing  program  must  "find"  the  appropriate  boundary 
conditions  that  satisfy  the  partial  differential  equation  (l)  and  the  surface 
energy  balances. 


Trajectory  Equations 


For  flight  we  must  solve  simultaneously  the  conservation  equations  and 
the  trajectory  equations  of  motion.  ¥e  use  the  two-dimensional  trajectory 
equations  (with  variable  mass)  for  entry  in  a  meridional  plane,  as  shown  below 
(ref.  18). 


(52b) 


(53a) 


(53b) 


dt  S]2 


(54) 


We  also  use  a  hypersonic  approximation  that  considers  the  ambient  atmospheric 
enthalpy  to  have  a  constant  value: 


(55a) 


Eaa  =  0.00836 


(55b) 


where  is  the  effective  average  constant  atmospheric  enthalpy  in  cal/g, 
and  E38  (an  input  to  the  computing  program)  has  the  units  km^/sec^.  For 
Earth  entries,  O.5  has  been  used  for  Ess. 
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Equations  (53)  and  (54)  can  "be  solved  numerically  in  a  straightforward 
manner.  As  programmed^  has  been  taken  as  constant  for  a  given  body. 

The  atmospheric  scale  height^  has  been  programmed  so  that  it  may  change  at 
selected  values  of  atmospheric  density,  D.  One  computing  option  for  Earth 
entries  provides  automatic  changes  in  through  three  values  to  represent 

the  AEDC  atmosphere.  For  arbitrary  planet  entries,  four  arbitrary  successive 
values  of  S-^  can  be  used.  (See  appendix  D.) 

The  quantity,  M/CpA  (g/cm^),  must  be  evaluated  for  use  in  equation  (53)- 
The  variation  of  M/A  has  been  set  up  empirically  as  a  function  of  the  sur¬ 
face  recession,  X.  This  is  equivalent  to  assuming  a  geometry  for  the 
recession  shape: 


M 

A 


(56) 


The  variation  of  the  drag  coefficient,  Cj),  through  the  free -molecule,  transi¬ 
tional,  and  continuum  regimes,  is  represented  by  a  bridging  equation  developed 
in  appenciix  C: 


where 


CpFM  -  Cpc 
Cpc 


(58) 


and  Eg  depends  on  the  body  shape.  The  free-molecule  drag  coefficient, 

CpFM^  niay  given  the  value  2.  The  parameter,  E14,  has  some  dependence  on 
body  shape  (and  flight  conditions),  but  will  often  be  assigned  the  value  zero 
(the  value  for  a  sphere  in  air).  We  combine  equations  (56)  and  (5T)  'to  obtain 

1  +  E13X  +  C3  ^e  ^ 

=  J - —  (59) 

CpA  VCpcVi 

In  this  equation,  Cjjq  is  shown  grouped  with  the  initial  value  of  M/Cj)qA, 
because  this  initial  quantity  is  used  in  the  computing  program.  The  empiri¬ 
cism  in  equation  (59)  can  also  be  used  to  account  for  any  change  in  Cpc  with 
change  of  body  shape. 


Miscellaneous  Relations 

Nose  radius.-  The  effective  nose  radius,  R,  is  calculated  as  a  quantity 
that  has  an  empirical  variation  with  the  surface  recession. 
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R  = 


1  +  EiyX^  +  Cl 


1 


(60) 


It  is  desirable  to  evaluate  the  constants  in  equation  (6o)  from  experimental 
data.  Otherwise^  one  must  estimate  the  geometrical  shape  of  the  ablation 
surfaces. 


Density  ratio  across  normal  shock. -  In  several  of  the  relations  presented^ 
the  density  ratio  across  a  normal  shocks  P2i^  appears.  For  wind-tunnel  calcu¬ 
lations^  P21  is  considered  to  remain  constant  for  a  particular  case.  The 
value  of  P21  is  a  portion  of  the  input  data  for  these  cases.  For  flight 
calculations^  P21  changes  and  must  be  continuously  calculated.  For  flight 
cases  we  use  the  equations: 


■Ei2RV^D 


2-V- 
12e  ^ 


1  +  0.08  log. 


0.1226' 


10  V  D 


Equations  (61)  and  (62)  are  empirical  relations  that  give  good  fits  to  data 
for  air.  Equation  (61)  accounts  for  nonequilibrium  effects^  while  equa¬ 
tion  (62)  fits  equilibrium  air  data^  such  as  presented  in  reference  I9.  The 
equations  provide  a  valid  approximation  for  any  gas  mixture  that  consists 
predominantly  of  nitrogen.  The  entire  analysis  is  not  overly  sensitive  to  the 
evaluation  of  p^^^.  It  appears  as  a  square  root  in  equations  (ll)  and  (35) 
which  are  approximations  themselves.  Values  of  En  and  E12  that  have  been 
used  for  Earth  entries  are  1.0  and  0.0001^  respectively. 


Other  Calculated  Quantities 

Although  not  always  required  in  the  analysis  described  in  the  preceding 
subsections^  several  other  quantities  are  calculated  because  they  are  of 
interest.  These  are  described  below. 


Thermal  thicknesses .  -  We  calculate  the  thermal  thickness^  based  on  the 
wall  temperature  gradient. 

This  is  meaningful  only  when  (8T/8y)^  <  0  (otherwise  A  is  arbitrarily 
assigned  the  value  10®) .  We  also  calculate  a  viscosity  thickness^  which 
we  define  as  the  depth  at  which  [x  =  ep^.  The  corresponding  viscosity  thick¬ 
ness  temperature  is  Ta,-  The  quantity,  Ta.»  is  determined  from  Tt^t  and  the 


and  the 


ness  temperature  is  The  quantity^  determined  from  Tw  and  the 

viscosity  representation  formula,  equation  {85) ,  belov.  The  quantity,  A|_i,  is 
determined  hy  interpolation  of  temperature  profile  data.  If  >  0  or  if 

T/^  exists  nowhere  in  the  temperature  profile,  is  assigned  the  value  zero. 
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stored  energy  comparison  vlth  exponential  temperature  profile.-  Using  a 
quantity,  cp,  we  compare  each  calculated  temperature  profile  with  an  exponen¬ 
tial  profile  having  the  same  (ST/3y)^.  The  quantity  cp  is  defined  as  the 
ratio  of  energy  stored  to  the  stored  energy  associated  with  the  corresponding 
exponential  temperature  profile  (with  form  similar  to  equation  (A1)  in 
appendix  A) .  Both  energies  are  calculated  as  constant  property  approxima¬ 
tions  .  For  the  area  under  the  exponential  profile  we  have: 

Aexp  =  A(Tw  -  To)  (l  -  +  Rb/f)  a(T^  -  T^) 


Then, 


Tn)dy. 


A(T^  -  To) 

Removal  of  melted  material  by  pressure  gradient  and  surface  shear. -  We 
calculate  a  term  that  compares  the  removal  of  melt  by  the  pressure  gradient 
and  the  surface  shear.  We  define  the  quantity,  an  approximate  evalua¬ 

tion  of  twice  the  ratio  of  the  portion  of  the  melting  velocity,  ux,  due  to  the 
pressure  gradient  to  that  due  to  surface  shear .  The  exact  quantity  would  be 
twice  the  ratio  of  the  first  term  to  the  second  term  in  equation  (9),  at 
y  =  0.  However,  we  use  the  approximation,  p  «  p^e  ^  h',  in  the  integrals  and 

obtain  the  simple  approximate  expression, 

2P"A, . 

a,d  .  --r^  (65) 


Surface  recession  due  to  vaporization  and  melting*  -  We  are  interested  in 
comparing  the  surface  recession  due  to  vaporization  vith  that  due  to  meiting* 
We  define  F  as  the  instantaneous  ratio  of  the  surface  recession  rate  due  to 
vaporization  to  the  total  surface  recession  rate  (due  to  vaporization  and 
melting) . 

F  =  3^  (66 

•Vbf 

For  the  ratio  of  the  surface  recession  due  to  vaporization  to  the  total  sur¬ 
face  recession  at  any  time,  t: 


Jlti _ 

pt 

J  hBFl'^’*' 


l-vwjdti 


_  ^vap 


Flow  lines  in  the  material.-  We  are  interested  in  the  path  of  the  flow  of 
material  up  to  the  point  of  melting  or  vaporizing.  This  is  of  particular 
Interest  in  the  study  of  tektites.  To  obtain  the  flow  lines,  we  make  an 
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approximate  quasi -steady  state  analysis. 

dx  =  ux  dt 

dy  =  V  dt 

dy/dx  =  v/ux 


Xr(y^t)  = 


(u/v) dy^ 


where  x^p  is  any  selected  small  value  of  x  at  y^^i  (where  u  =  G) .  Then 
X  =  XpX;gp  gives  x(y)  at  time  t  for  a  flow  line. 

Aerodynamic  deceleration.-  For  flight  cases^  we  are  interested  in  the 
aerodynamic  deceleration  of  the  hody.  We  calculate  the  ratio  of  the  aero¬ 
dynamic  force  to  the  mass  of  the  body  to  obtain  the  aerodynamic  deceleration. 
We  define  a^g  as  the  component  of  deceleration  due  to  aerodynamic  drag^ 
normalized  by  the  gravitational  acceleration  of  the  planet: 


^Dg  *■ 


M  \  /gp 

VlO^. 


The  absolute  value  of  the  total  normalized  aerodynamic  acceleration  is 


1  +  (L/Dr) 


Reradiation  and  apparent  emissivity.-  The  terms  for  the  fluxes  of 
radiation  at  the  front  and  back  surfaces  are  given  in  equation  (76) .  The  rate 
of  reradiation^  from  the  front  surface  (absolute  value)  is  Fj^3.  For  opaq.ue 
bodies^ 

%s  ==  (71) 

The  expression  for  the  reradiation  from  the  back  surface  of  opaque  bodies  has 
a  similar  form.  For  opaque  material^  the  surface  emissivities^  and 
will  be  assigned  values.  For  semitransparent  material  representation^  the 
evaluation  of  is  as  given  in  equation  (45) • 

For  transparent  bodies^  the  situation  is  somewhat  different.  Here^  we 
wish  to  calculate  effective  surface  emissivities  (which  will  vary  with  time). 
The  flux  of  reradiation  from  the  front  surface  is  given  by 

=  'Ir  -  (72) 

Then  we  calculate: 
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(73) 


'FF  -  oT/ 


Similarly,  for  the  hack  face; 


F(yBF^’t) 


^BF  aTT3^4 


The  quantities,  e^rp  and  €gpi,  in  equations  (73)  (7^)  fh.®  calculated 

effective  surface  emissivities  for  transparent  material. 

Energy  Balance 

Equation  (l)  is  solved  by  numerical  finite  difference  methods,  and  it  is 
desirable  to  check  the  accuracy  of  solutions  obtained.  This ^ is  done  by  cal¬ 
culating  a  group  of  energy -integral  terms  listed  belo*w,  summing  them  up,  and 
determining  the  residual  (error)  in  the  sum.  The  magnitude  of  the  energy 
integrals  is  also  of  interest,  because  it  shows  the  disposition  of  the 
energies  involved;  this  knowledge  gives  an  insight  into  the  processes  of 
ablation. 

We  make  the  listing  of  the  energy  rate  terms  as  follows: 

The  total  convective  heat -transfer  rate  into  the  material  is 

^con  ^  BE 

The  net  radiative  heating  rate  into  the  material  is 


=  F(0,t)  -  F(y^^,t)  +  (1  -  Bi6)[(1|^  -  “  "^o  ^  ^ 


where 


Bi6  =  1  for  the  transparent  case 

Bis  =  0  for  the  opaque  and  semitransparent  cases 


F  =  0  for  the  opaque  and  semitransparent  cases 

The  energy  accounted  for  by  the  rate  of  vaporization  (positive  into  the  mate 
rial)  is 
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The  energies  put  into  the  material  are  accounted  for  hy  terms  that  involve 
convection  of  the  ablating  material  in  the  x  and  y  directions  and  hy  a 
storage  term.  For  the  rate  of  increase  of  stored  energy^ 


^stor  ^  dt 


yBF  pT(yi) 


Cp  dTg  dy^ 


The  X  direction  convection  energy  rate  term  (positive  out)  vd.ll  be 


%.con 


r^BF  - 

/  u  /  Cp  dTs  dy 

Jn  J’V 


The  y  direction  convection  energy  rate  term  (positive  out)  is 


'-'rp 

■^o 

The  error  in  the  energy  rate  balance  will  be  a  residual  term^  ^res* 

Sres  ~  ‘^con  ‘^rad  'Ivap  "  ‘Istor  ”  '^ucon  "  ‘^vcon 

The  residual^  ‘Ires^  shovrs  the  accuracy  of  the  energy  rate  balance  at  any  time 

t. 

The  cumulative  energy  balance  shows  the  total  size  of  the  various  terms 
involved  and  the  error  accumulation.  We  compute  the  following  integrals: 


q  dti 
^con 


'Irad 


(82b) 


^irnTi  <iti 


(82c) 


'stor  P 


ypF^^)  pT(y  ^t)  pY^p(t±)  pT'(y  .tj^) 

/  Cp  dTs  dy  -  p  /  /  Cp  dTs  dy 


(82d) 
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'^con  “  ‘lucon 


(82e) 


^Ivcon 


(82f) 


So,  for  the  accumulated  residual  we  have: 


Qres  '^con  '^rad  *^vap  "  '^stor  “  *^ucon  "  *^vcon 


METHOD  OF  THE  HOMERICAL  PROGRAM 


Representations  of  Physical  Properties 

Most  of  the  equations  in  the  previous  section  contain  (at  least 
implicitly)  quantities  representing  physical  properties  of  the  ablating  mate¬ 
rial  or  external  gas.  The  density  of  the  ablating  material,  p,  is  assumed  to 
be  constant,  but  the  other  pertinent  physical  properties  are  considered  to  be 
temperature  dependent.  The  temperature  dependence  of  the  physical  properties 
has  been  left  unspecified;  however,  to  obtain  numerical  solutions,  it  will  be 
specified  by  formulas  with  constants  that  can  be  arbitrarily  chosen  and  read 
into  the  computing  program.  The  representations  used  are: 


Equilibrium  vapor  pressure  of  ablating  material,  atm 


I>ve  =  ® 


t!+^3 


Viscosity  of  ablating  material,  poises 


p  =  e 


B4 

T-B14 


Specific  heat  of  ablating  material,  cal/g 


c  =  B0  +  EiT  - 


Thermal  conductivity  of  ablating  material,  cal/cm  sec 


K  =  Bg  + 


+  E20T' 


Average  specific  heat  of  external  gas,  cal/g  (see  ec^.  (l4)) 


Cp  =  Eio  +  EsT-^ 
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Transformation  of  Coordinates  and  Finite  Differencing 

As  described  in  the  ANALYSIS  AM)  METHOD  OF  SOLUTION  section^,  the 
procedure  in  solving  the  system  of  equations  presented  is  to  make^  in  effect^ 
all  equations  auxiliary  to  equation  (l)  which  is  solved  with  its  appropriate 
boundary  conditions  as  given  in  equations  (46)  and  (47)  or  (51)*  Equation  (l), 
a  partial  differential  equation,  is  solved  by  a  finite  difference  scheme;  some 
complications  are  produced  by  the  decrease  of  length  with  time  between  the 
front  and  back  surfaces  due  to  front  surface  recession.  In  this  work,  it  was 
elected  to  keep  the  number  of  grid  points  constant,  so  the  distance  between 
grid  points  was  allowed  to  shrink  with  decreasing  length  of  the  ablator.  This 
was  accomplished  by  the  following  transformation  of  independent  variables  from 
t,  y  to  s,  Ti,  where 

s  =  t  (89) 


_ 

^  L  -  X(t) 

Then, 

^  _ 

Sy 


=  A_  +  (dX/ds)Ti  5 

Bt  Ss  L  -  x(s)  3:1 


(92) 


This  transformation  alters  somewhat  the  form  of  equation  (l)  and  the  other 
equations  as  actually  put  into  the  numerical  computing  program. 

Equation  (l)  is  solved  numerically  by  an  explicit  (forward  difference) 
scheme.  In  finite  difference  form,  the  partial  derivatives  of  the  tempera¬ 
ture,  T,  are  represented  as: 


(93) 


I 


2k 


o 


n+i 

n 


o 


s 


o  o  o 


m-l 


Sketch  (h) 


where  m  -  m,  m  +  1  are  grid  point 
numbers  on  the  ti  (depth)  scale  and  n, 
n  +  1  are  numbers  on  the  s  (time) 
scale  as  shown  in  sketch  (b) .  Finite 
increments  of  s  and  t]  are  indicated 
by  the  A  symbol. 


Stability  and  Accuracy  of  the  Finite  Difference  Equation 


In  solving  a  parabolic  partial  differential  equation  by  a  forward 
01‘ence  scheme^  there  is  always  a  stability  requirement  to  be  met.  For  the 
finite  differencing  of  the  transformed  version  of  equation  (l),  the  stability 
requirement  turns  out  to  be: 


Z  = 


As 

L 

^  fK\ 

At 

/  K  \ 

L(Ati)2J 

L  -  X(s)_ 

ypc J  " 

_(Ay)2_ 

(96) 


The  stability  parameter,  Z,  is  printed  out  by  the  computing  program  for  each 
grid  point  at  each  time  printed.  The  increments  of  Aq  remain  constant  with 
time,  but  the  increments  of  Ay  decrease  with  time,  as  indicated  by  equa¬ 
tion  (90) ^  until  ablation  is  concluded.  Thus,  Z  tends  to  increase  somewhat 
while  ablation  is  proceeding,  and  this  must  be  considered  in  selecting  Initial 
increments.  Since  Z  should  not  exceed  l/2,  it  is  not  possible  to  use  the 
computing  program  to  calculate  to  the  point  of  complete  extinction  of  an 
ablating  material.  When  the  present  numerical  program  is  used  for  transparent 
material,  the  storage  limitation  requires  that  the  initial  value  of  the  ratio 
of  length  (depth)  to  the  smallest  At]  be  <  I665  (see  spacing  sketch  in 
appendix  D ) . 

A  gross  check  on  the  accuracy  of  numerical  solutions  obtained  is  provided 
by  the  running  energy  balance  and  the  cumulative  energy  balance  (see  MALYSIS 
AND  METHOD  OF  SOLUTION  section) .  An  additional  check,  standard  in  numerical 
work,  can  be  made  by  varying  As  and  At)  and  noting  the  resultant  variations 
produced  in  the  numerical  solutions.  This  check,  in  effect,  determines  the 
adequacy  of  representing  the  derivatives  by  difference  quotients  with  the 
finite  increments  as  chosen  (-eqe.  (93)^  (9^)j  (95))* 


Boundary  Conditions  for  Transparent  Material 

It  was  noted  in  the  ANALYSIS  AND  METHOD  OF  SOLUTION  section  that  the 
computing  program  can  be  arranged  to  modify  the  surface  energy  balance  equa¬ 
tions  (46)  and  (47)  for  the  transparent  case.  Equations  (46)  and  (4t)  are 
rigorous  as  written,  but  they  require  a  fine  spacing  of  grid  points,  and 
therefore  small  time  increments  for  stability.  Particularly  near  the  front 
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face^  the  derivative  of  the  radiation  flux^  can  change  hy  orders  of 
magnitude  in  a  short  distance.  With  a  "normal”  spacing  of  grid  points^  the 
fine  structure  of  the  variation  of  F  and  g  are  not  well  represented^  and 
this  can  result  in  energy  balances  that  are  not  accurate  (large  residuals). 
In  the  computing  program^  the  front  surface  energy  balance  that  ve  actually 
use  for  the  transparent  case  (Bis  =  l)  is: 


pb-vVv  +  ”  '^1/2) 


(97) 


where  F1/2  is  the  radiation  flux  midway  between  the  front  surface  and  the 
first  interior  finite -differenced  grid  point  (point  number  1  +  K2  in  FORTRAN 
terminology;  see  spacing  sketch  in  appendix  D) .  For  the  back  surface  energy 
balance, 


57)3^  ■  %r  -  '^37(Pbf  -  (98) 

where  ^BF-i/^  radiation  flux  midway  between  the  back  surface  grid 

point  and  the  nearest  interior  grid  point.  The  quantities,  E36  and  E37,  are 
constants  read  into  the  program  and  will  have  values  >  0  and  <  1.  When 
E36  =  E37  =  0,  we  recover  equations  (46)  and  (4?)  exactly.  E36  and  E3Y  can 
be  adjusted  to  give  optimum  energy  balances;  very  good  energy  balances  have 
been  obtained  with  E36  =  E37  =1. 


ILLUSTRATIVE  EXAMPLES 


Examples  shown  below  illustrate  the  use  of  the  numerical  computing 
program  as  applied  to  several  types  of  ablation.  Calculated  and  measured 
results  are  compared  for  all  examples  except  the  last  one,  which  is  a  Martian 
entry.  The  disposition  of  energies  for  typical  examples  given  is  summarized 
in  table  I. 


Tektite  Glass  in  a  Wind  Tunnel 

Tektite  glasses,  ablated  at  high  enthalpies  in  an  arc -jet  wind  tunnel, 
furnish  examples  of  ablators  that  both  vaporize  and  melt.  Typical  comparisons 
between  calculated  and  measured  values  of  surface  recession  and  surface 
brightness  temperature  are  shown  in  figures  1  and  2.  The  measured  points  in 
figure  2  are  actually  a  spread;  measurements  on  other  tektite  glasses  fell 
between  these  points .  The  agreement  can  be  seen  to  be  very  good,  which  lends 
confidence  that  data  for  flights  involving  tektite  glass  can  be  successfully 
calculated  (ref.  l) .  In  both  figures  the  calculations  were  made  for  a  trans¬ 
parent  glass  and  for  a  semitransparent  glass,  and  there  is  little  difference 
between  the  results  of  the  two  methods  of  computation.  The  glasses  used  in 
these  examples  ablate  by  melting  more  than  by  vaporizing  because  of  the 
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Figure  1. -  Comparison  of  calculated  and  measured 
surface  recession  of  tektite  glass  in  a  wind 
tunnel. 


t,  sec 


Figure  2.-  Comparison  of  calculated  and  measured 
brightness  temperature  of  tektite  glass  in  a 
wind  tunnel. 


moderate  enthalpy  in  the  wind  tunnel,  the  low  viscosities,  and  low  vapor 
pressures  of  the  glasses.  For  the  glasses  in  figure  1,  ahout  1  percent  of  the 
ablation  is  due  to  vaporization,  and  for  the  glasses  in  figure  2,  vaporization 
accounts  for  less  than  1  percent.  At  higher  enthalpies  the  relative  amount  of 
vaporization  Increases .  Table  I  gives  the  disposition  of  energies  calculated 
for  the  semitransparent  glass  of  figure  1.  Because  of  the  small  amount  of 
vaporization,  very  little  heat  is  blocked,  and  most  of  the  Incoming  energy  is 
accounted  for  by  melting. 


Tektite  Entry  Calculation 


The  results  of  calculating  an  entry  for  a  typical  opaque  tektite  are 
shown  in  figures  3  and  k.  An  entry  speed  of  11.0  km/sec  and  an  entry  angle  of 
-30°  were  used  for  the  calculations.  These  conditions  correspond  to  a  typical 
deduced  trajectory  for  a  Victoria  australlte  (ref.  2,  fig.  22).  Figure  3 
shows  the  calculated  values  of  velocity,  surface  temperature,  surface  reces¬ 
sion,  and  surface  recession  due  to  vaporization.  The  free -molecule  and  con¬ 
tinuum  regimes  are  also  distinguished.  Time  zero  is  arbitrarily  selected  as 
"far  out"  before  any  appreciable  aerodynamic  heating  has  begun.  This  example 
illustrates  the  response  of  a  material  that  vaporizes  readily,  with  about 
Ik  percent  of  the  ablation  due  to  vaporization  and  the  rest  to  melting.  As 


Va,i  =  ii.o  km /sec 


0  .2  .4  .6  .8  1.0  1.2  1.4  1.6 

y,  cm 


Figure  3.  Calculated  variation  of  surface  tern-  Figure  4.  -  Calculated  temperature  profiles  for  a 
perature,  velocity^  total  ablation,  and  Victoria  australite  entering  Earth’s  atmo- 

vaporized  ablation  for  a  Victoria  australite  sphere;  Rj_  =  O.816  cm. 

entering  the  Earth’s  atmosphere;  Rj_  =  O.816  cm. 
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the  tektite  heats  up^  its  surface  begins  to  ablate  at  a  temperature  in  excess 
of  2000  The  surface  temperature  and  ablation  rate  reach  a  maximum  and 

then  fall  off  as  the  body  slows  down.  The  end  of  ablation  occurs  rather 
abruptly^  and  the  remainder  of  the  flight  is  that  of  a  solid  body  being  aero- 
dynamically  cooled.  Measurements  of  the  amount  of  ablation  at  the  stagnation 
point  (refs.  1  and  2)  on  recovered  tektites  yield  values  not  greatly  different 
from  the  7.3  mm  calculated  for  this  example.  The  calculations  indicate  that 
for  this  flight,  a  negligible  portion  of  the  ablation  occurred  in  the  free- 
molecule  regime.  The  portion  of  the  ablation  in  the  transitional  regime  was 
2k  percent,  compared  to  the  majority  of  ablation  in  the  continuum  regime 
(76  percent) .  For  smaller  tektites  and  shallow  entry  angles  the  percent  of 
ablation  in  the  free-molecule  and  transitional  regimes  will  be  greater;  for 
large  vehicles  this  portion  of  ablation  is  generally  small. 

In  figure  k  are  shown  variations  of  the  calculated  temperature  profiles 
and  surface  recession  at  selected  values  of  time.  This  figure  gives  a  fairly 
complete  picture  of  the  internal  heating  and  eventual  cooling  of  the  body 
during  its  flight.  The  rise  in  the  back  temperatures  is  due  to  base  heating. 
Measurements  on  recovered  bodies  by  photoelasticity  techniques  show  locked -in 
thermal  stresses  that  vary  from  a  depth  of  0.2  to  0.35  cm,  corresponding  to 
the  calculated  depth  of  0.23  cm  (ref.  2').  The  post-flight  observations  of 
thermal  stresses  and  deduced  surface  recession  on  recovered  bodies  are  com¬ 
patible  with  the  calculated  results.  The  energy  disposition  for  the  flight 
calculated  is  given  in  table  I.  The  vaporization  that  occurs  in  this  flight 
causes  substantial  heat  blockage;  the  bulk  of  the  energy  is  accounted  for  by 
heat  blockage,  heating  and  vaporizing,  and  heating  and  melt  flow,  these  three 
quantities  being  of  about  the  same  order  of  magnitude. 


Reentry  Flight  With  Silica  Glass  Heat  Shield 


A  calculation  was  made  for  a  reentry  flight  of  a  nose -cone  with  an  opaque 
silica  glass  heat  shield.  The  vehicle  and  flight  are  described  by  Hidalgo  and 

Kadanoff  in  reference  20.  For  this  trajec¬ 
tory  and  this  heat -shield  material  about 
Ik  percent  of  the  ablation  was  due  to  vapori¬ 
zation  which  is  comparable  to  the  tektite 
entry  case  previously  discussed.  The 
recovered  reentry  vehicle  allowed  the  amount 
of  ablation,  X,  at  the  stagnation  point  to  be 
determined.  The  physical -property  inputs  in 
this  case  correspond  to  opaque  silica 
(ref.  20),  but  both  the  transparent  and  the 
semitransparent  options  of  the  computing 
program  were  run  with  the  results 


X( transparent)  2.  11 

X(observed) 


X( semitransparent) 


X(observed) 


=  1.15 
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The  corresponding  ratio  as  calculated  by  Hidalgo  and  Kadanoff  using  their 
quasi-steady  ablation  analysis  was  about  1.10  and  by  Chapman  and  Larson 
(ref.  l) ^  using  an  integral  method  of  calculation^  0.92.  This  illustrates 
that  the  amount  of  ablation  on  simple  materials  such  as  glass  can  be  computed 
to  the  order  of  10-  to  15-percent  accuracy  by  several  methods,  including  two 
of  the  options  of  the  present  computing  program.  In  view  of  the  inevitable 
angle -of -attack  variations  in  flight,  which  cause  the  stagnation  point  of 
maximum  heating  to  wander  somewhat  over  the  nose  and  thus  reduce  somewhat  the 
maximum  recession,  the  observed  difference  between  calculated  and  measured 
ablation  is  in  the  expected  direction.  The  energy  disposition  for  both  the 
transparent  and  semitransparent  calculations  is  shown  in  table  I.  It  is  of 
interest  that  the  two  calculations  yield  energy  proportions  that  are  nearly 
the  same,  although  the  internal  temperature  distributions  are  different 
because  of  radiative  transmission  in  the  transparent  case.  In  both  cases  the 
total  ablation  is  moderate,  so  that  the  actual  amount  of  vaporization  is 
moderate  and  the  heat  blockage  term  is  relatively  small. 


Teflon  Model  in  an  Arc-Jet  Wind  Tunnel 


Under  normal  ablative  conditions,  tetraf luoroethylene  polymer  (Teflon) 
undergoes  a  surface  depolymerization  and  vaporization  of  the  monomer  at  a  sur¬ 
face  temperature  of  approximately  7^0  ®K.  There  is  no  one  specific  tempera¬ 
ture  at  which  the  reaction  occurs,  but  a  sharply  rising  reaction  rate  with 
temperature  in  this  neighborhood  essentially  controls  the  surface  temperature 
of  an  ablating  model  (refs.  21,  22).  Under  these  conditions,  the  viscosity  of 
Teflon  remains  high,  so  the  process  can  be  said  to  resemble  a  sublimation 
(with  the  reaction  rate  determined  by  an  Arrhenius  type  rate  equation) .  In 
performing  the  calculations  for  Teflon  ablation,  it  was  assumed  that  any 
energy  involved  in  possible  chemical  reactions  between  the  Teflon  vapor  and 
the  external  gases  could  be  neglected. 


Comparisons  between  calculated  and  experimentally  measured  surface 
recession  for  Teflon  are  shown  in  figure  5.  The  experiments  by  G.  Lee  and 
R.  Sundell  (ref.  23)  were  performed  in  an  arc- jet  wind  tunnel  for  four  values 
of  enthalpy.  It  is  seen  that  the  agreement  obtained  is  quite  good,  with  the 
possible  exception  of  the  7^0  Btu/lb  total  enthalpy  case.  It  is  thought  that 


Figure  5-"  Comparison  of  calculated  and  measured 
surface  recession  of  Teflon  in  wind  tunnel. 


the  generally  satisfactory  agreement 
shown  in  the  figure  indicates  the 
validity  of  the  method  calculation  and 
also  that  the  physical  properties  of 
the  substance  have  been  adec^uately 
represented.  The  front  face  mass  loss 
rate  was  essentially  reaction  rate  con¬ 
trolled  for  the  four  cases  shown  in  the 
figure . 

The  disposition  of  the  calculated 
energies  for  the  2000  Btu/lb  total 
enthalpy  case  is  shown  in  table  I. 

The  heat  blockage  term  is  fairly  large 
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because  all  the  ablated  material  leaves  in  the  vapor  state.  The  largest  term 
is  the  energy  for  heating^  depolymerizing^  and  vaporizing. 


Teflon  Heat  Shield  in  a  Mars  Entry 


The  calculations  for  this  example  illustrate  the  application  of  the 
computing  program  to  an  entry  vith  a  proposed  Mars  probe  (ref.  24).  A  spheri¬ 
cal  capsule  of  6l.0-cm  diameter  has  been  assumed^  with  a  1-cm  thick  Teflon 
heat  shield^  entering  the  Martian  atmosphere  in  an  oriented  attitude.  Four 
hypothetical  atmospheres  were  assumed  as  tabulated  below. 


Composition  (vol.) 


Scale  height)  km 


100^  Ns 

7.8 

91%  Ns,  %  COs 

7.8 

100^  Ns 

20.0 

91%  Ns,  9%o  COs 

20.0 

Subsequent  to  making  these  calculations^  data  obtained  from  the  19^5  Mariner 
occultation  experiment  have  indicated  that  the  scale  height  of  the  Martian 
atmosphere  is  about  9  km  (ref.  25).  The  assumed  atmospheres  with  the  7-8  km 
scale  heights  thus  appear  to  be  the  more  realistic  ones.  Calculations  were 
made  for  an  entry  velocity  of  7.92  km/sec  and  for  two  entry  angles,  -90^  and 
-20*^.  An  M/Cj)A  for  continuum  flow  of  3*9^  g/cm^  was  assumed  for  the  vehicle, 
and  m/a  was  held  constant  while  the  Cp  varied  through  the  transition  from 
the  free -molecule  to  the  continuum  regimes. 

The  -90^  entries  have  the  greater  peak  heating  rates,  but  the  -20*^ 
entries  absorb  more  total  heat  and  are  the  more  critical  from  a  heat -shield 
standpoint.  The  heat-shield  responses  are  compared  using  calculated  values 
of  the  stagnation-point  recession  as  tabulated  below. 


Atmosphere 


Entry  angle, 
deg 


Total  stagnation 
point  recession^  cm 


N2  small-scale  height 
H2  large-scale  height 
N2-CO2  small-scale  height 
N2-CO2  large-scale  height 
N2-CO2  small-scale  height 
^2~C02  large-scale  height 


-90 

0.098 

-90 

.159 

-90 

.121 

-90 

.179 

-20 

.196 

-20 

•  3i^l 

Of  the  -90^  entries  in  the  table,  the  most  severe  environment  would  be  the 
mixed  atmosphere  and  the  large-scale  height  (although  the  small-scale  height 
appears  to  be  more  realistic).  For  a  given  scale  height,  the  mixed  atmosphere 
gives  somewhat  more  ablation  than  the  nitrogen  atmosphere  because  there  is 
more  radiation  from  the  carbon  dioxide  (ref.  I7) .  The  time  of  exposure  to 
heating  is  roughly  proportional  to  scale  height  for  two  otherwise  similar  tra¬ 
jectories,  and  the  total  heat  absorbed  is  approximately  proportional  to  the 
square  root  of  exposure  time  (and  therefore  scale  height).  This  approximate 
relationship  between  scale  height  and  total  recession  for  a  given  atmospheric 
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composition  can  be  deduced  from  the  table.  The  strong  dependence  of  the 
heat-shield  response  on  any  uncertainty  of  knowledge  of  atmospheric  scale 
height  is  of  interest,  since  this  trend  will  presumably  apply  to  any  heat- 
shield  material  and  any  planet. 


The  disposition  of  energy  for  the  -20°  entry  with  the  mixed  atmosphere 
and  short -scale  height  is  shown  in  table  I.  As  with  the  wind-tunnel  results 
for  Teflon,  the  two  large  energy  terms  are  the  heat  blockage  term  and  the 
term  that  accounts  for  heating,  depolymer izing,  and  vaporizing.  In  the  envi¬ 
ronment  of  the  -20°  Martian  entry,  the  considerable  amount  of  ablation  of 
material  to  the  vapor  state  accounts  for  the  very  large  heat  blockage  term. 


MARS  ENTRY  ASSUMED  ATMOSPHERE 

ENTRY  VELOCITY  =  7.92  km/sec  91%  Ng,  9%  CO2  (VOL) 

ENTRY  ANGLE  =*-20“  SCALE  HEIGHT  =  20.0  km 

TOTAL  MATERIAL  LOSS  =  1.56  kg  =  3.43  lb 


Figure  6.-’  Approximate  calculation  of  response 
of  a  Teflon  heat  shield  in  a  Mars  entry. 


As  an.  illustration^  the  analysis 
has  also  been  used  in  an  approximate 
manner  to  calculate  the  quantities  of 
interest  around  the  front  hemisphere  of 
the  spherical  capsule  for  the  -20'^ 
entry  in  the  mixed  atmosphere  with  the 
(probably  overly  severe)  large  scale 
height.  The  results  are  summarized  in 
figure  6  which  shows  the  variation  of 
total  recession  and  front  and  back  face 
maximum  temperatures  around  the  hemi¬ 
spheric  heat  shield  as  well  as  the 
total  mass  loss  for  this  hypothetical 
case.  These  results  illustrate  how 
variable  material  thickness  may  be  used 
in  heat -shield  design. 


TABLE  I.-  TYPICAL  ENERGY  BALANCES  (percentages) 


Tektite  glass 
(semitransparent) 
wind  tunnel 
(fig.  1) 

Tektite 
Earth  entry 
(opaque) 

(figs.  3,4) 

Silica  glass 
heat -shield 
Earth  entry 
(transparent) 

Silica  glass 
heat -shield 
Earth  entry 
(semitransparent) 

Teflon  wind 
tunnel  (fig,  5) 
hs=2000  Btu/lb 
Pg=0.34  atm 

Teflon 

Mars  entry 
7i=-20° 

Vcol=7.92  km/seo 

Convection  In 
(hot  wall) 

100.0 

100.0 

100.0 

100.0 

100.0 

100.0 

Heat  "blocked 

2.9 

31.9 

11.4 

11.5 

35.4 

71.9 

(Net  convection  in) 

(97.1) 

(68.1) 

(88.6) 

(88.5) 

(6lt.6) 

(28.1) 

Net  radiation  in 

-13-3 

-10.8 

-28.4 

-27.6 

-4.8 

2.4 

Heating  and 
vaporizing 

5.2 

24.9 

21.5 

21.2 

50.1 

25.5 

Heating  and  melting 

69.1 

27.2 

28.6 

31.6 

0 

0 

Stored 

4.6 

6.3 

10,4 

6.9 

9.9 

5.1 

Error 

4.9 

-1.1 

-0.3 

1.2 

-0.2 

-0.1 
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CONCLUDim  REMAEKB 


A  generalized  analysis  of  stagnation-point  ablation  has  been  presented 
for  solving  a  variety  of  problems  involving  melting  and  vaporizing^  subliming^ 
or  surface  chemical  reactions.  The  flexibility  of  the  analysis  has  been 
demonstrated  through  the  presentation  of  several  varied  illustrative  examples. 
In  general^  it  is  expected  that  accuracy  of  answers  obtained  will  depend 
largely  on  the  degree  of  knowledge  of  the  physical^  chemical^  and  thermody¬ 
namic  properties  of  the  ablating  material^  as  these  are  necessary  inputs  for 
the  computing  program.  The  procedure  of  relating  calculations  for  a  given 
material  to  experiment  wherever  possible  lends  confidence  to  calculations  for 
the  same  material  exposed  to  other  conditions  which  cannot  be  verified  by 
observation. 


Ames  Research  Center 

National  Aeronautics  and  Space  Administration 
Moffett  Fields  Calif July  22^  I966 

129-03-12-01-00 
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APPENDIX  A 


PRINCIPAL  NOMENCLATURE 


In  performing  computing  machine  calculations^  some  purely  FORTRAN 
q^uantities  are  used^  particularly  among  input  data^  which  have  no  counterpart 
among  the  synbols  listed  below.  These  quantities  are  in  appendix  D^  wherein 
all  FORTRAN  quantities  are  listed. 

a  y  direction  body  force  per  unit  mass^  cm/sec^  (acceleration  in 

flight) 

ar^  component  of  deceleration  due  to  aerodynamic  drag  normalized  by 

^  gravitational  acceleration  of  planet^  dimensionless  (eq.  (69)) 

ag  absolute  value  of  aerodynamic  acceleration  normalized  by  gravita¬ 

tional  acceleration  of  planet^  dimensionless  (eq.  (70)) 


^cq^^cv^ 

Ar»m 


frontal  area^  cm^ 

free -molecule  acconmodation  coefficients  for  heat  transfer^  mass 
loss^  X  momentum  (for  surface  shear respectively^  dimension¬ 
less 

corrected  mass-loss  accommodation  coefficient  (eq.  (39c) 
dimensionless 


constant,  defined  by  equation  (27) 

constant,  defined  by  equation  (ll) 

constant,  defined  by  equation  (32a) 

constant,  defined  by  equation  (15) 

melt-off  parameter,  defined  by  equation  (65) ^  dimensionless 
Sutherland  constant,  (eq.  (C38)) 

Arrhenius  frequency  factor^  cm/sec  (eq.  (39c)) 
constant  in  vapor  pressure  (eq.  (84)) 
constant  in  vapor  pressure  (eq.  (84)) 
constant  in  viscosity  (eq.  (85)) 
constant  in  viscosity  (eq.  (85)) 
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Be 

By 

Bg 

Bii 

Bi4 

Bi6 

C 

c" 

Cp 

Cp 

CppM  / 
Cl 
02 
C3 
C4 
Cs 
D 

®aZ 

®rr, 

Eip 

E2(z) 

El 

E2 


constant  in  specific  heat  (eq.  (86)) 

constant  in  specific  heat  (eq.  (86)) 

constant  in  thermal  conductivity  (eq,  (87)) 

constant  in  convective  heat  blockage  factor  (eqs.  (28)^  (30)) 

(See  discussion  following  eq.  (C7),) 

constant  in  viscosity  (eq.  (85)) 

constant;  Big  =  1.0  for  transparent  case;  Bie  =  0  for  opaque  and 
semitransparent  cases  (eqs.  (46),  (47),  (51)) 

specific  heat  of  body  material,  cal/g 

heat  capacity  per  unit  area  of  backing  material^  cal/cm^ 

(eq.  (51)) 

specific  heat  of  a  gas  at  constant  pressure,  cal/g 

average  specific  heat,  external  gas,  cal/g  defined  by  equation  (l4) 

drag  coefficient,  continuum  drag  coefficient,  free-molecule  drag 
coefficient,  respectively,  dimensionless 

constant  in  nose  radius  (eq.  (60)) 

constant  in  nose  radius  (eq.  (60)) 

constant  in  M/A  (eq.  (56)) 

constant  in  M/A  (eq.  (56)) 

constant,  vorticity  correction  in  equation  (15) 
free -stream  density,  g/m^ 

allowable  error  in  (selected),  ^K;  allowable  disagreement 
between  ^w  obtained  from  equations  (l)  and  (46) 

error  in  after  last  iteration,  *^K;  will  be  <  e^^^ 

Arrhenius  activation  temperature,  (eq.  (39^)) 

exponential  integral  (second  degree),  defined  in  equation  (43) 

constant  in  specific  heat  (eq.  (86)) 

constant  in  thermal  conducitivity  (eq.  (87)) 
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Ea 

E4 

Es 

Eg 

E7 

Eg 

E9 

Eio 

Ell 

Ei2 

Ei3 

Ei4 

Ei5 

Ei6 

Ei7 

Eis 

Ei9 

E20 

E33 

E35 

^36^^37 

E38 


constant  In  average  specific  heat  (eq.  (88) ) 
constant  in  gas-cap  radiation  (eq.  {hk)) 
constant  in  gas-cap  radiation  (eq.  (44)) 
constant  in  gas-cap  radiation  (eq.  (44)) 

constant  to  account  for  shift  of  vaporization  equilihrium 
(eq.  (26)) 

constant  in  expression  for  shear  blowing  factor  (eq.  (34)) 

constant,  defined  by  equation  (58) 

constant  in  average  specific  heat  (eq.  (88)) 

constant  in  expression  for  (eq.  (61)) 

constant  in  expression  for  (eq.  (61)) 

constant,  defined  by  equation  (50) 

constant  depending  on  body  shape  and  flight  conditions  in  drag 
bridging  (eq.  (57)) >  also  equation  (C46) 

constant  accounting  for  body  force  in  wind  tunnel  expression  for 
P"  (eq.  (lOa)) 

constant  used  in  tumbling  correction  (eq.  (l2)) 
constant  in  nose  radius  (eq.  (60)) 
constant  in  M/A  (eq.  (56)) 

constant,  defined  by  equation  (49)  as  the  ratio  of  turbulent  to 
laminar  base  heating 

constant  in  thermal  conductivity  (eq.  (87)) 

constant  used  in  expression  for  front  face  emlssivity  for  semi¬ 
transparent  body  (eq.  (45)) 

constant,  asymptotic  value  of  \lr  (eq.  (28) ) 

(See  discussion  following  eq.  (07)*) 

constants  used  in  equations  (97)  and  (98)  to  modify  surface  energy 
balances  for  the  transparent  case 

average  ambient  enthalpy  for  flight  case,  km^/sec^  (eq.  (55)) 
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F 

^RS 

'^1/2. 

^BF-1/2 

F 

h 

g 

h 

hs 

llY 

h 

h 

K 

Ka 

Ka 

K 

Ktu 

L 


ratio  of  actual  viscosity  to  undissociated  (Sutherland)  value^ 
dimensionless 

radiation  flux,  cal/caP  sec  (eq.  (4l)) 

reradlation  rate  from  front  surface,  cal/cm^  sec  (eqs.  (7I),  (72)) 

radiation  flux  midway  "between  front  surface  and  first  interior 
finite -differenced  grig  point,  point  number  1  +  k2  in  FORTRAN 
terminology;  see  Spacing  Sketch  in  appendix  D;  cal/cm^  sec 
(eq.  (97)) 

radiation  flux  midway  Between  Back  surface  grid  point  and  nearest 
interior  grid  pointy  cal/cm^  sec  (eq.  (98)) 

ratio  of  surface  recession  rate  due  to  vaporization  to  total  sur¬ 
face  recession  rate,  dimensionless  (eq.  (66)) 

ratio  of  surface  recession  due  to  vaporization  to  total  surface 
recession,  dimensionless  (eq.  (67)) 

gradient  of  radiative  flux,  3F/6y,  cal/cm^  sec 

gravitational  acceleration  of  planet,  cm/sec^  (eq.  (53B)) 

enthalpy,  cal/g 

stagnation  enthalpy,  cal/g 

latent  heat  of  vaporization,  cal/g 

average  ambient  enthalpy  for  flight  case,  cal/g  (eq.  (55)) 

h/hg,  dimensionless 

(1  -  h)/(l  -  h^),  dimensionless 

thermal  conductivity,  cal/cm  sec 

mean  free  path  constant,  moles  (eq.  (c4o)) 

Reynolds  analogy  factor  (eq.  (31))^  dimensionless 

fractional  time  center  point  is  exposed  to  stagnation  conditions, 
dimensionless  (eq.  (l2)) 

correction  factor  for  heat  transfer  and  other  quantities  due  to 
oscillation  in  flight,  dimensionless  (eq.  (l2)) 


Le 
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initial  depth  of  material,  cm 
Lewis  number,  dimensionless 


L/Dr 

m 

» 

m 

M 

00 

M 

n 

P 

Pts 

Vy 

^ve 

^•vm 

P 

P’* 

a 

%c 

^FM 

vort=o 


lift/drag  ratio^  dimensionless 

molecular  weight 

mass  loss  rate^  g/cm^  sec 

mass  of  hody^  g;  also  grid  point  number  in  finite  difference 
computation 

free-stream  Mach  number^  dimensionless 
me/my^  dimensionless  (eg.  (-29)) 

index  of  refraction,  dimensionless;  also  exponent  in  equation  (30) 

pressure,  dynes/cm^  or  atm,  as  specified 

pressure  downstream  of  normal  shock,  atm 

actual  vapor  pressure,  atm 

equilibrium  vapor  pressure,  atm 

modified  equilibrium  vapor  pressure,  atm  (eq.  (26)) 

ratio  of  pressure  of  external  gas  to  modified  equilibrium  vapor 
pressure  of  ablated  vapor,  dimensionless  (eq.  (26)) 

negative  of  second  derivative  with  respect  to  x  of  external 

pressure  with  a  correction  for  body  force,  dynes/cm"^;  defined  in 
equation  (8);  evaluated  in  equation  (lO) 

heat-transfer  rate,  cal/cm^  sec;  also  dynamic  pressure,  dynes/cm^ 

surface  convective  (continuum)  heat -transfer  rate  with  no  blowing, 
cal/cm^  sec  (eq.  (l5)) 

surface  convective  (free  molecule)  heat -transfer  rate,  cal/cm^  sec 

(eq.  (17)) 

surface  convective  heat -transfer  rate  with  no  blowing,  bridged 
between  qQ^  and  qp]y[,  cal/cm^  sec  (eq.  (20)) 

qoo  corrected  for  tumbling  or  oscillation,  cal/cm^  sec  (eq.  (2l)) 

q^^  without  vorticity,  cal/cm^  sec;  C©  =  0  in  equation  (15) 

gas-cap  radiation  rate,  cal/cm^  sec  (eq.  (4^4-)) 
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a  surface  convective  (continuum.)  heat-transfer  rate  with  "blowing, 

cal/cm^  sec  (eq.  (l6)) 

q.  surface  convective  heat -transfer  rate  (all  regimes),  "bridged 

^  "between  q^^  and  q-pj^j,  cal/cm^  sec  (eq.  (l8)  or  (22h)) 

a  q^„  corrected  for  tumbling  or  oscillation,  cal/cm^  sec 

^  ^eq.  (19)  or  (23)) 

q  p.  initial  combined  convective  and  radiative  heating  rate, 

cal/cm^  sec  (eq.  (B3a)) 

Energy  transfer  rates  listed  below  pertain  to  an  energy  balance,  and, 
where  applicable,  are  combined  rates  for  front  and  back  surfaces. 

q  total  convective  heat- transfer  rate  into  material,  cal/cm^  sec 

(eq.  (75)) 

net  radiative  heating  rate  into  material,  cal/cm^  sec  (eq.  (76)) 

g  energy  rate  due  to  vaporization  (positive  if  energy  is  released 

into  the  material),  cal/cm^  sec  (eq.  (77)) 

g  .  rate  of  increase  of  stored  energy  in  the  material,  cal/cm^  sec 

(eg.  (78)) 

q  X  direction  convection  energy  rate  of  the  material  (positive 

out),  cal/cm^  sec  (eq.  (79)) 

g  y  direction  convection  energy  rate  of  the  material  (positive 

out),  cal/cm^  sec  (eg.  (80)) 

g^gg  residual  in  energy  rate  balance,  cal/cm^  sec  (eg.  (8I)) 

■s 

^con^^rad,  time  integrals  of  corresponding  q  values,  cal/cm^  (eq.  (82)); 

Qvap.?Qstor^  '  terms  in  total  energy  ‘balauce 

^con^^vcon 

residual  in  total  energy  balance,  cal/cm^  (eq.  (83)) 

R  nose  radius,  cm 

Rl^  body  radius,  cm  (eqs.  (48),  (50)) 

Rg ratio  of  base  to  front-face  laminar  convective  heating 

normalized  by  respective  enthalpy  potentials,  for  exposed 
back  surface,  dimensionless  (eq.  (48)) 


Reynolds  number 
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effective  coefficient  of  reflection  for  planar  radiation, 
dimensionless  (eq..  (4-2)) 

Reg  Reynolds  nuniber  based  on  enthalpy  velocity  (eq.  (C20)) 

R  universal  gas  constant,  ergs/mole  °K 

o 

Rp  planet  radius,  km  (eq.  (53)) 

s  transformed  time  coordinate,  sec  (eq.  (89)) 

S  collision  cross-section  area,  cm^  (eq.  (C4o)) 

atmospheric  scale  height,  km  (eq.  (5^)) 

initial  scale  height  for  entry  into  arbitrary  atmosphere,  km 
^  (eq.  (b6) ) ,  and  contained  in  equation  (BIO) 

8h  ^^h  successive  scale  heights  in  arbitrary  atmosphere,  km: 


Sh  = 

vhen 

Poo  —  P0O2 

Sh  =  Sh2 

vhen 

P0O2  ^  ^00 

Sh  =  Sh3 

when 

Poo  ^  P003 

time^  sec 

initial  time^  sec;  time 

at  which  front 

temperature,  °K 

brightness  temperature  (emissivity  unity) ,  °K 
reference  temperature,  °K 

viscosity  thickness  temperature,  °K;  temperature  at  which 
p  =  e.p^  (at  depth  Ap) 

horizontal  component  of  trajectory  velocity,  km/sec 

X  direction  velocity  of  external  gas  at  edge  of  boundary  layer, 
km/sec  (eq.  (ll)) 

velocity  of  material  in  x  direction,  cm/sec 
Up. 10®,  cm/sec 


,  sec 
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u* 


V 

V 


X 


Xp 

y 

y 

z 


a 

a2 

7 

6 

&X 

6* 

A 


€ 


A 

M- 


u/ue;  dimensionless  (eq.  (C2l) ) 

vertical  component  of  trajectory  velocity  (positive  upward)^  km/sec 
velocity  of  material  in  y  direction^  cm/sec 
surface  recession  rate^  cm/sec  (eg.  (6a)) 

enthalpy  velocity^  km/sec;  defined  as  =  O.OO836  hs  (eg.  (13)) 

free -stream  (trajectory)  velocity^  km/sec  (eq.  (52a)) 

longitudinal  coordinate  along  meridian^  cm 

flo¥  line  ratio,  dimensionless;  defined  "by  eguation  (68) 

transverse  coordinate  normal  to  surface  (inward),  cm 

boundary-layer  transverse  coordinate,  cm;  appendix  C 

stability  parameter  for  finite  differencing,  dimensionless; 

(eg.  (96));  must  be  <  l/2 

absorption  coefficient,  internal  radiation,  cm"^  (eg.  (Al)) 
absorption  coefficient,  gas-cap  radiation,  cm"^  (eg.  (Al)) 
trajectory  angle,  deg,  positive  above  horizontal  (eg.  (52b)) 
increment  (appendix  C) 

mass  fraction,  dimensionless  (e-g.  (CIO)) 
displacement  thickness,  cm 

thermal  thickness,  cm  (eg.  (63));  also  increment  e.g., 

2I1  =  enthalpy  potential,  cal/g 

unspecified  characteristic  boundary- layer  thicknesses,  cm 
(appendix  C) 

viscosity  thickness,  cm;  depth  at  which  [x  =  e.p^ 

surface  emissivity,  opague  or  semitransparent;  effective  emissivity, 
transparent  (egs.  (73) ^  (7^))^  dimensionless 

transformed  y  coordinate,  cm  (eg.  (90));  also  dummy  variable 

mean  free  path,  cm 

viscosity,  poise 


ho 


f  Pr 
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density  of  ablating  material  (constant),  g/cm^ 

density  ratio  across  normal  shock,  dimensionless j  for  wind  tunnel 
cases,  assigned;  for  flight  cases ,  equations  (6l)  and  (62) 

free-stream  atmospheric  density  in  Earth  sea-level  atmospheres, 
dimensionless,  D/1226 

values  assigned  to  at  which  changes  of  scale  height  in  an 

arbitrary  atmosphere  occur;  see  ^hs 

Stefan  constant,  1.369x10"^^  cal/cm^  sec  °K^;  also  Prandtl  number, 
dimensionless 

shear,  dynes/cm^ 

ft 

approximation  of  the  ratio  of  stored  energy  to  the  stored  energy 
associated  with  an  exponential  temperature  profile,  dimensionless 
(eq.  (6i|-)) 

surface  recession^  cm 

surface  recession  due  to  vaporization,  cm;  Xya,p  =  (^T)) 

characteristic  recession  depth,  cm,  used  in  tumbling  correction 

(eq.  (12)) 

convective  heat  blockage  factor,  dimensionless  (eqs.  (I6),  (28)) 
modified  convective  heat  blockage  factor,  dimensionless  (eqs.  (22), 

(24)) 


Subscripts 

a  actual 

BP  back  face 

c  continuum 

cha  change  of  wind  tunnel  conditions 

d  diffusion 

e  external  gas  or  outer  edge  of  boundary  layer 

eq  equilibrium 

FF  front  face 
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FM 


free  molecule 


i  initial;  this  subscript  can  be  combined  with  the  others 

max  maximum 

o  no  blowing 

oc  no  blowing^  continuum 

off  shut  off  of  wind  tunnel 

r  reverse 

ref  reference 

s  stagnation  (or  settling  chamber) 

u  undissociated 

V  vapor  expelled 

w  wall  (front  face) 

wc  wall^  continuum 

wd  wall^  diffusion 

wFM  wall^j  free  molecule 

1  dummy  var  i  ab  le 

2  behind  normal  shock;  also  average  condition  between  shock  wave  and 

body  (appendix  C) 

free  stream 

Superscript 

’  X  derivative 
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APPENDIX  B 


EQUATIONS  FOR  STARTING  VALUES 


To  start  the  solution  to  equation  (l),  it  is  necessary  to  assign  initial 
conditions.  These  will  normally  consist  of  a  relatively  low  temperature  pro¬ 
file  which  can  exist  before  the  onset  of  ablation.  The  particular  selection 
of  initial  conditions  is  generally  not  critical  as  their  influence  damps  out 
in  a  short  time.  The  initial  temperature  profile  that  we  assume  is  an  expo- 
nential  type. 


Ti  (y)  =  To  +  (Ti 


-y/Ai  ^ 


-(yBF-y)/^i 


If  T  •  is  selected  near  Tq,  this  profile  amoTonts  to  a  small  perturbation 
on  the^constant  Tq  profile.  We  can  take  the  y  derivative  of  Tj^  at 
y  =  0,  equate  it  to  the  ratio  of  the  heat  flux  (eq.  (46))  to  thermal  conduc¬ 
tivity  at  the  wall^  and  solve  for  the  initial  thermal  thickness^  Ai- 


"  G.vi  +  "  Bi6)qRi 


(B2a) 


We  use  the  more  simple  approximation: 


Kwi(Twi  -  ^q) 

q„i  +  (1  -  BisjqRi 


(B2b) 


where  Bis  =  1  transparent  case  and  Bis  =  0  Bhe  opaque  and  semi¬ 

transparent  cases.  As  described  below  in  the  section^  Flight  Cases ^  one 
option  allows  Aj_  to  be  assigned  a  value  instead  of  obtaining  it  from  equa¬ 
tion  (B2b).  The  initial  convective  and  radiative  heat  fluxes  needed  in  equa¬ 
tion  (B2b)  are  obtained  differently  for  wind  tunnel  or  flight  as  described^ 
below.  In  calculating  q^  and  qj^^,  the  initial  free -stream  <iensitp  Di,  is 
needed.  For  wind  tunnel  calculations,  Di  will  be  known;  for  flight,  D^  can 
be  assigned  or  calculated  as  described  below. 


Ablation  in  a  Wind  Tunnel 


The  starting  of  a  wind  tunnel  is  visualized  as  a  sudden  step  of  a  Beat^ 
flux.  We  define  a  combined  initial  heat  flux,  q^Ri^  ‘’^Be  sum  of  the  Initia 
convective  and  radiative  fluxes. 


‘^wRi  ”  '^wi  ‘^Ri 


(B3a) 


^3 


We  assume  that  t  =  1  and  that  the  convective  heat  transfer  is  in  the 
continuum  regime  calculated).  With  =  ^oc±  have: 

^wRi  ~  loci  ^  ^hi  (^3'b) 

We  evaluate  from  equation  (15)  and  from  equation  (44)^  with 

""  using  I'yi  for  li*  I'wi  small^  loci  S'^id  therefore  Iv/Ri 

are  approximately  constant  for  a  short  period  of  time  (eq.  (l5))-  The  classi¬ 
cal  conduction  problem  with  a  constant  heat  flux  and  constant  properties 
(ref.  26 y  p.  56)  can  be  used  as  an  approximation  for  this  case  to  determine 
the  time  at  which  the  front  face  temperature  arrives  at  the  assumed  T^i. 

This  turns  out  to  be: 


ti  = 


(Tyi  "  Tq) 

^IwRi^ 


(b4) 


One  can  set  ^  greater  value  is  not  necessary^  but  it  gives  the 

computing  program  a  smooth  start. 


Flight  With  Arbitrary  Initial  Conditions 

In  starting  the  flight  calculations^  we  will  use  an  assigned  initial 
velocity  V^j_  and  a  flight  path  angle ^  at  (arbitrary)  time^  t^  =  0.  The 

initial  atmospheric  density^  I)±  (equivalent  to  a  starting  altitude)^  can  be 
assigned^  as  well  as  an  initial  thermal  thickness^  Aj^.  With  an  assumed  ^y±y 
the  initial  profile  is  determined  from  equation  (B1) .  Using  this  starting 
procedure^  we  do  not  require  that  Aj_  be  consistent  with  the  relationship 
given  in  equation  (B2b) .  However^  we  calculate  q.^j_^  q^^y  and  q^j^j_  as  fig¬ 
ures  of  interest  as  described  below  in  the  section^  Initial  Conditions  for 
Entry  Flight . 


Initial  Conditions  for  Entry  Flight 

An  alternative  starting  procedure^  valid  for  an  entry  flighty  is  to  use 
an  assigned  and  and  an  assumed.  T^q^  and  to  calculate  the  entry  into 

an  exponential  atmosphere  [l  =  which  "will  raise  the  front  face 

temperature  to  the  assumed  T^i.  The  convective  heating  during  this  initial 
^  part  of  the  entry  will  be  considered  to  be  of  the  free -molecule  type  (and  we 
do  not  calculate  Qoci^’  write  (eq.  (l9))* 

^wi  " 

where 

^tvi  ^FMi 
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The  evaluation  of  in  equation  (l2)  specializes  to 


Ktu  =  K 


so  we  have 


We  evaluate  from  equation  (17)^  and  using  equations  (B3a)  and  (44)^  we 

have 


KAc^DiVi 


‘IwRi  =  CpT^i)  +  KE4RiDi^^Vi^® 


(B5) 


where  we  use  as  an  approximation  for 


The  quantity  q-^j^i  is  the  heating  rate  at  time  t;j^  =  0.  Up  to  time 
ti  =  0,  the  heat  flux  will  he  approximately  an  exponential  function  of  time^ 
through  the  exponential  variation  of  D.  This  is  similar  to  a  classical  prob¬ 
lem  in  heat  conduction  (ref.  26^  p.  45)  which  yields  an  exponential  tempera¬ 
ture  profile.  In  equation  (B5)  we  can  replace  hy  D  as  an  exponential 

function  of  altitude;  we  can  approximate  the  slightly  varying  enthalpy  poten¬ 
tial  as  a  constant  (the  value  in  eq.  (B5));  and  we  can  integrate  the  heating 
flux  over  time  from  t  =  -00  to  t  =  tj_  =  0.  We  obtain^  then^  the  total  heat 
absorbed: 


Q 


ShiKA,cqDi(Vi^  -  0.00836  CpT^i)  ShiKE4RiDi 


E5  E0-1 


total 


0.0836  I  sin  7^ 


(b6) 


Es  Sin  7. 


We  can  also  approximate  the  total  heat  absorbed  as 


When  we  substitute  the  profile  equation  (B1)  into  equation  (B7); 


^total  ^  P^wi  ('^wi  ”  TQ)Aj_(l  +  1  -  e 


■yBpAi^ 


(B7) 

we  obtain 

(B8a) 


or,  approximated  further, 

^total  ^  P^wi(Tvi  -  (B8b) 

We  substitute  q^^.  and  q„.  (the  first  and  second  terms,  respectively,  on  the 
right  side  of  eq.  (B5)J  into  equation  (B2b)  for  Aj^,  put  this  into  equa¬ 
tion  (B8b),  and  have  finally 


4 


ir\ 


Q 


total  ~ 


(B9) 


E5,.  Ee 


0.0836"  '  °'°°^36  CpT^i)  +  KE4RiDi  % 


We  now  eliminate 
tion  of  the  form 


*^total 


from  equations  (b6)  and  (B9)  and  we  obtain  an  equa- 


t  KiDi 


(BIO) 


where  Ki  and  K2  are  constants  and  is  the  only  unknown. 

The  procedure,  then,  for  starting  the  entry  calculations  is  as  follows. 

We  assume  a  and  we  know  the  scale  height,  Sh^,  far  out  in  an  atmosphere. 

We  use  equation  (BIO)  to  calculate  Dp;  we  can  then  calculate  and  as 

the  first  and  second  terms,  respectively,  on  the  right  side  of  equation  (B5)- 
Finally,  we  obtain  Aj_  from  equation  (B2b),  and  we  put  into  the  profile 

equation  (B1) .  We  can  see  that  the  assumption  of  T^p  fixes  the  Dp,  or,  in 
effect,  fixes  the  altitude  at  which  we  start  time  zero.  For  this  case, 

T^p  >  Tq  is  necessary  in  order  to  have  a' finite  starting  altitude. 
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APPENDIX  C 


BRIDGING  BETWEEN  FREE -MOLECULE  AND  CONTINUUM  REGIMES 


For  a  number  of  applications  or  situations,  it  is  known  that  the 
transitional  regime  Between  free-molecule  and  continuum  flow  must  he  consid¬ 
ered  The  Bridging  formulas  used  in  this  work  for  the  transitional  regime  are 
presented  without  derivations,  in  the  ANALYSIS  ^  METHOD  OF  SOLUTION  section 
L  equations  (l8),  (20),  {2k),  (37),  (i+O),  and  (57)  •  The  derivations  of  these 
equations,  Based  principally  on  a  simple  kinetic  theory  model,  are  given  in 
this  appendix.  These  Bridging  formulas  replace  previously  used  equations 
which  were  of  purely  empirical  form  (ref.  l) . 

Front  Face  Normal  Velocity 

For  the  front  face  normal  velocity  or  mass  loss  rate,  a  Bridging 
relationship  is  required  Between  the  free-molecule  and  continuum  regimes.  We 
will  consider  first  the  free-molecule  or  reaction  rate -controlled  regime,  and 
we  will  distinguish  Between  the  evaporation  or  suBlimation  case  and  the 
chemical  reaction  case. 

For  the  case  involving  evaporation  or  sublimation,  ve  write  the  Langmuir 
equation  (eq.  (76)  of  ref.  27)  for  the  mass  loss  rate  into  a  vacuum  as; 

%M  =  -^cvP^e^dy j  2rtRgT^ 

where  the  constant,  C^jy,  is  the  pressure  of  a  standard  atmosphere  in 
dynes/cm^  so  that  Pyg"^  is  measured  in  atmospheres,  my  is  the  molecular 
weight  in  g,  and  Rg  is  the  universal  gas  constant  in  erg/mole  °K.  Equa¬ 
tion  (Cl)  is  Based  on  the  rate  of  molecules  crossing  a  unit  area,  in  this 
case  Impacting  against  a  unit  surface  area,  and  the  accommodation  coefficient, 
Acv,  is  the  fraction  of  the  molecules  that  stick  to  the  surface  and  condense. 
At  equilihrium  the  evaporation  and  condensation  rates  are  equal  so  there  is 
then  no  net  evaporation  rate.  At  a  given  temperature,  the  rate  of  surface 
impacts,  and  therefore  the  rate  of  condensation,  is  proportional  to  the  actual 
vapor  pressure  aBove  the  surface,  Py.  So  we  can  write  for  the  condensation 
rate: 

Sr  =  ^ 


and  for  the  net  rate  of  evaporation  (or  suBlimation) 


m  — 


17 


(C3b) 


If  the  equilibrium  vapor  pressure  is  shifted  or  modified  by  the  presence  of 
other  gaseous  materials  in  the  boundary  layer^  then  the  largest  actual  vapor 
pressure  at  the  surface^  p^^  that  can  be  reached  is  the  modified  equilibrium 
vapor  pressure,  Py^^i  (see  eq.  (26)),  rather  than  p^e-  accordingly  modify 
equation  (C3b)  and  use  the  approximate  expression: 


(C3c) 


For  the  chemical  reaction  case  we  rewrite  the  Arrhenius  rate  equation 
for  a  first  order  reaction  (eq.  (39©))  as 


mpM  =  pSe 


-Et/Tw 


When  the  modified  equilibrium  vapor  pressure,  V-Ym.?  exists  above  the  surface, 
the  net  reaction  rate  is  zero,  or  the  reverse  reaction  rate  equals  the  forward 
reaction  rate.  At  a  given  temperature,  the  rate  of  impact  of  vapor  molecules 
with  the  surface  will  again  be  proportional  to  the  actual  vapor  pressure  at 
the  surface,  p^.  We  can  assume  that  the  reverse  reaction  rate,  is  propor¬ 

tional  to  the  rate  of  surface  impacts  and  therefore  to  the  actual  vapor  pres¬ 
sure,  Py,  and  we  again  have  equation  (C3c)  for  the  net  reaction  rate. 

We  now  consider  the  continuum  or  diffusion  controlled  regime  for  the 
front  face  mass  loss  rate.  The  limiting  value  of  the  diffusion  controlled 
front  face  normal  velocity,  jv^cU  given  in  equation  (38),  and  we  define 
%  =  p|v^c|  (obtained  by  putting  p  on  the  left  side  of  eq.  (38)),  as  the 
limiting  value  or  maximum  mass  diffusion  rate.  As  noted  in  the  discussion 
following  equation  (38),  this  evaluation  is  obtained  using  P  and  ^  which 
ultimately  depend  on  the  modified  equilibrium  vapor  pressure  at  the  surface, 
p^.  As  distinct  from  the  theoretical  or  limiting  value,  we  next  consider 
the  actual  mass  transfer  rate  by  diffusion.  The  actual  rate  can  be  shown  to 
be  approximately  proportional  to  the  actual  vapor  pressure  at  the  surface, 
p^,  by  inserting  actual,  rather  than  limiting,  evaluations  of  P  and  \|/  into 
equation  (38)  .  In  these  actual  evaluations,  P  will  be  based  on  p^  rather 
than  Py^  and  the  \|/  asymptote  is  assumed  to  be  0  with  Bn  pa  1.  We  then 
obtain 

(C5a) 


This  is  straightforward  when  <  P-^  •  When  p-y^  >  Ptg^  3.  modification  is 

necessary  because  in  equation  (38)  we  naTe,  in  effect^  given  p^  an  upper 

limit  of  0.999^999  Ptg* 


(C5i3) 


Let 


^  /pW\ 


(c6) 


and  we  have  equation  (G5a)  with  replacing  1%. 

If  a  quasi-steady  state  is  assumed;,  there  is  no  accumulation  or  depletion 
of  vapor  in  the  boundary  layer.  Thus  the  actual  mass  loss  rate  calculated  hy 
the  free-molecule  (reaction-rate)  method^  and  the  actual  mass  loss  rate  calcu¬ 
lated  as  a  diffusion  process  must  he  the  same;  physically  this  means  that  ^ 
the  actual  vapor  pressure  at  the  surface,  p.^.,  must  have  a  value  such  that  m 
obtained  from  equations  (C3c)  and  (C5)  will  be  the  same.  We  eliminate  Pv/Pvm 
between  equations  (C3c)  and  (C5a)  and  obtain 


1 

S 


%M  % 


(C7) 


At  high  mass  loss  rates,  one  can  surmise  that  diffusion  control  may  merge  into 
a  hydrodynamic  control  with  an  interf  ace  between  the  two  fluids .  In  this  case, 
equation  (38)  may  yield  values  of  the  "diffusion"  rate  that  are  too  large  at 
the  high  mass  loss  rates.  Finite  values  of  the  "diffusion"  rate  will  be 
obtained  from  equation  (38)  when  the  asymptote  E35  =  0,  and  Bn  =  1. 

For  the  (somewhat  unusual)  situation  when  p^jj^/p^j^g  >0.8,  it  is  suggested 
that  one  use  E35  =  0  and  Bn  isd  1;  this  insures  that  =  pI'^'wcU  calcu¬ 
lated  using  equation  (38),  does  not  become  unrealistically  large.^  For  the 
condition  with  p-^m  >  we  should  actually  replace  with  in 

equation  (C7).  However,  this  situation  wiU  generally  be  one  of  rate  control 
in  which  %]y[  will  be  small  relative  to  ma.  (whose  calculated  value  may  be 
too  large)  and  will  be  still  larger.  As  an  approximation  for  all 

conditions,  we  will  accordingly,  use  equation  (C7)  with  m<j. 

When  we  cancel  out  the  constant  density,  p,  from  equation  (CT),  we  have 
equation  (hO)  for  the  front  face  normal  velocity,  Vy.  A  somewhat  similar  line 
of  reasoning  for  the  evaporation  case  is  in  reference  28  although  bridging 
equations  are  not  presented.  These  bridging  forms  (eqs.  {hO)  and  (C7))  are 
considered  to  be  valid  approximations  over  the  complete  spectrum  from  rate^ 
control  to  diffusion  control;  the  use  of  the  bridging  relationship  automati¬ 
cally  places  control  in  the  proper  regime.  Some  verification  of  this  bridging 
relationship  is  given  in  figure  5  hy  the  comparisons  between  calculations  and 
experiment  for  the  ablation  of  Teflon  in  a  wind  tunnel.  These  examples  are 
essentially  reaction  rate  controlled,  except  for  hg  =  3OOO  Btu/lb  which  is 
considered  to  be  in  the  transitional  regime. 
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Convective  Heat  Transfer 


To  calculate  heat  transfer  in  the  transitional  regime^  we  require  a 
bridging  between  the  continuum  and  free -molecule  heat  transfer  values  given  by- 
equations  (l6)  and  (17).  The  bridged  result  is  shown  in  equation  (18).  A 
very  simple  first-collision  model  is  used  in  the  analysis.  A  typical  packet 
of  free  molecules  is  assumed  to  enter  the  boundary  layer  and  make  a  first 
collision  with  molecules  already  there.  The  free  molecules  then  become  part 
of  the  continuum  boundary  layer  with  average  energy  equal  to  that  of  the 
boundary  layer  at  the  point  of  the  collision.  The  energy  given  up  by  the  free 
molecules  on  collision  is  assumed  to  be  ultimately  taken  up  by  the  wall  (by 
successive  collisions  in  the  continuum  boundary  layer  and  impact  with  the 
wall),  since  there  is  no  piling  up  of  energy  in  a  quasi-steady  state  boundary 
layer.  (Some  of  the  free  molecules  will  make  their  first  collision  with  the 
wall  and  give  up  energy  directly.) 


We  use  a  normalized  enthalpy,  h  = 
in  sketch  (d).  We  define  an  effective 


Sketch  (d) 


entering  the  boundary  layer  and  making 
allows  us  to  write: 


h/hg,  and  the  coordinate  system  shown 
collision  thickness,  such  that 
collisions  occurring  within  have 

an  effect  on  wall  heat  transfer,  while 
collisions  occurring  outside  of 
have  a  negligible  effect  on  wall  heat 
transfer;  is  thus  a  kind  of  bound¬ 
ary  layer  thickness.  We  can  rewrite 
equation  (I7)  as 

qpM  =  Ki(DV^)hs(l  -  EJ  (C8) 

where  Ki  contains  the  conversion  of 
units  and  the  accommodation  coefficient, 
Aeq  (assumed  to  be  approximately  unity). 
Then,  for  a  small  packet  of  molecules 
its  first  collision,  our  assumed  model 


[lisKi(DVj]5X(l  -  h)  (C9) 

where  6X  is  the  mass  fraction  of  molecules  that  make  the  first  collision 
between  Y  and  Y  +  dY  (or  between  E  and  h  +  dh)  in  the  boundary  layer;  bqp]y[ 
is  heat  given  up  by  the  fraction  6X.  In  using  Ki  in  equation  (C9) ,  we  are 
assuming  that  the  fraction  of  energy  given  up  in  the  molecular  collisions  is 
the  same  as  the  fraction  given  up  by  collisions  with  the  wall  (approximately 
unity) ;  this  is  thought  to  be  within  the  framework  of  approximations  being 
made  in  the  derivation.  The  mass  fraction,  6X,  is  evaluated  in  terms  of  the 
mean  free  path,  A  (ref.  29^  eq.  (103-7) )• 


6X 


-Y/A  ^ 
®  A 


(CIO) 
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The  mean  free  path.  A,  is  actually  a  function  of  Y  since  A  varies  inversely 
with  density  in  the  boundary  layer.  We  will  consider  A  to  have  some  con¬ 
stant  mean  value  in  the  hoimdary  layer  in  order  to  perform  an  integration 
(below) .  This  assumption  appears  to  be  within  the  framework  of  approximations 
being  made.  In  accord  with  our  quasi- steady-state  assumption,  we  write 


all  8X 


and  using  equation  (C9)  we  have 


(Cll) 


all  SX 


[hsKi(DV^) ]  [hsKi(DV^) ] 


=  ^  8X(l  -  h) 


all  8X 


(C12) 


We  insert  5X  from  equation  (CIO)  and  have: 


[hoKi(DVj] 


all  6X 


(C13) 


We  -will  use 


1  -  h...  VAi. 


(C14) 


and  ve  can  write: 


[hsKi(DVj] 


=  (1  -  hv) 


Ap -YAffi.,  -AiA 


(015) 


The  last  term  in  the  bracket  accounts  for  heat  transfer  from  the  molecules 
whose  free  path  is  greater  than  Aq  (ref.  29,  eq.  (IO3-8)).  Using  equa¬ 
tions  (C8  and  Cli)-)  we  have  finally 


(C16) 


We  will  use  the  simplest  form  for  f(Y/Aq)  in  equation  (Cl4)  (see  sketch  (d)). 


(C17) 
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We  perform  the  integration  in  equation  (Cl6)  and  obtain 


1  -  e 


-Z^/A 


(C18) 


To  consider  a  variation  in  the  ratio ^  visualize  a  change  in 

ambient  density^  The  mean  free  path^  A^  varies  as  and  Z\^  (being  a 

kind  of  boundary  layer  thickness)  varies  as  so  the  ratio^  variei 

as  In  the  free-molecule  limits  A/Aq^  becomes  very  large^  and  in  the 

continuum  limit,  very  small.  When  A/Ag  is  large,  q^^  approaches 
equation  (Cl8),  as  it  should.  When  A/^  is  small,  q^^  must  approach 
Using  equation  (Cl8)  ve  can  write  for  small  A/Z\q 


Z\q 


(C19a) 


(C19b) 


The  quantities,  q^j^  and  q^^^  are  calculated  values  for  given  (the  same)  condi¬ 
tions.  Since  the  characteristic  thickness,  Z\q,  has  not  been  specified,  it  can 
now  be  given  the  value  that  satisfies  equation  (C19b)  for  the  conditions 
imposed.  The  mean  free  path.  A,  has  been  considered  to  be  some  mean  value  in 
the  boundary  layer,  and  it  can  be  assigned  a  convenient  value,  say  As^  with 
Ag  required  to  satisfy  equation  (C19b)  with  this  A.  It  can  be  shown  that 
the  form  of  equation  (C19b)  is  a  consistent  relationship  by  making  use  of 
equations  (5*5)^  (6.23)^  (6.25)^  (6.^6)  of  reference  30^  relating  A  to 
a  characteristic  thickness  (displacement  thickness,  6*). 

We  substitute  equation  (C19b)  into  (Cl8)  and  we  have  equation  (l8)  as 
the  bridging  equation  for  convection  heat  transfer. 


Equation  (l8)  has  the  form  of  a  mo  no  tonic  function  of  and  it  has 

the  correct  free-molecule  and  continuum  regime  asymptotes.  The  derivation 
has  been  essentially  performed  from  the  free-molecule  end,  and  the  agreement 
with  the  continuum  end  has  been  forced.  In  the  derivation,  functional  forms 
for  f(y/Ag)  other  than  the  simple  one  chosen  (eq.  (CI7))  can  be  used.  The 
result  will  be  similarly  behaved,  but  more  complicated,  bridging  equations. 
It  is  possible  that  the  more  complicated  equations  that  can  be  obtained  may 
give  better  agreement  with  measured  data  for  some  specific  types  of  heat 
transfer . 


The  alternate  forms  of  the  bridging  relations  (eqs.  (20)  to  (24))  are 
obtained  directly  from  equation  (18)  as  outlined  in  the  MALYSIS  AM)  METHOD  OF 
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SOLUTION  section.  The  alternate  forms  yield  the  q.uantities^  and  ^ 

(which  are  mainly  of  conceptual  interest) . 

Comparisons  of  calculated  bridging  with  experiment  are  available  for  the 
case  with  no  blowing.  Experimental  data^  obtained  from  figure  5  of  refer¬ 
ence  31^  are  shown  plotted  in  figure  7  and  compared  with  calculated  values  of 

continuum  heat  transfer^  free -molecule 
heat  transfer^  and  the  bridged  heat 
transfer  value  obtained  from  equa¬ 
tion  (20) .  The  data  shown  are 
stagnation -point  heat -transfer  measure¬ 
ments  on  a  spherical  body  at  nominal 
Mach  numbers  of  5.7  and  8  with  stagna¬ 
tion  temperatures  of  2100^  to  2300^  R. 
The  quantities^  q^^  (eq.  (l5))^  IpM 
(eq.  (17))^  and  q^^  (eq.  (20))  are  all 
normalized  using  the  continuum  heat- 
transfer  rate  with  zero  vorticity^ 

Figure  7.-  Comparison  of  calculated  convective  loc  V  l^lis  normalizing  factor  iS 
heat  transfer  "bridging  with  experiment  for  VOrt=0 

the  case  without  blowing.  Obtained  from  equation  (l5)  with  the 

•vorticity  correction  dropped  (Cs  =  O) .  The  Reynolds  number  used  as  abscissa, 
Reg,  is  based  on  the  enthalpy  velocity  (eq.  (iSb))  and  the  stagnation  gas 
properties,  as  originally  used  in  reference  3I. 


Reg  - 


Vp.iDR 

lOpg 


(C20) 


In  the  experiments  of  reference  3I,  the  Reynolds  number  was  varied  by  mainly 
varying  the  free -stream  density,  D,  and  the  nose  radius,  R.  It  is  seen  from 
the  figure  that  the  calculated  bridging  checks  well  with  experiment. 

Comparisons  were  also  made  (not  shown  here)  with  measured  data  from 
reference  32  for  subsonic  heat  transfer  from  spheres.  Equation  (20)  checks 
these  data  reasonably  well. 

It  is  concluded  from  the  comparisons  made  that  the  bridging  relation  in 
equation  (18)  should  be  a  useful  approximation  for  the  general  case  with  blow¬ 
ing  since  equation  (20)  is  simply  a  specialization  (with  the  form  unchanged) 
of  equation  (18),  and  equation  (18)  has  the  correct  asymptotes  for  the  general 
case.  Experimental  verification  for  the  blowing  case  would  be  desirable. 


Surface  Shear 

For  the  calculation  of  the  (x  derivative  of)  surface  shear  in  the 
transitional  regime,  we  bridge  between  equations  (35) and  (36)  to  obtain  equa¬ 
tion  (37).  The  first  collision  model  used  is  the  same  as  that  used  for  con¬ 
vective  heat  transfer,  except  that  we  are  now  concerned  with  the  transfer  of 
X  momentum  rather  than  energy. 
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Sketch  (e) 


¥e  use  the  coordinate  system  shown 
in  sketch  (e)^  and  we  introduce 


(dUe/dx)^^Q 

We  assume  that  collisions  occurring 
within  an  effective  collision  thick¬ 
ness^  (a  kind  of  boundary  layer 
thickness)^  affect  surface  shear ^  while 
those  occurring  outside  have  a 

negligible  effect.  Making  use  of  equa¬ 
tion  (ll);  we  rewrite  equation  (36)  as 


(C21) 


(C22) 


where  K4  contains  the  conversion  of  units  and  the  accommodation  coefficient, 
Acm^  considered  to  be  approximately  unity.  Using  the  same  reasoning  and 
analogous  development  as  that  used  for  convective  heat  transfer,  we  can  write 

Ti-  ^  (C23a) 

all  5X 


^  SX(l  -  u*) 


(C231) 


all  BX 

We  evaluate  5X  according  to  equation  (CIO)  and  obtain: 


T*  =  ri 


V  'wFM 
We  use  (see  sketch  (e)) 


^A. 


U  Uq 


(1  -u.)e-^Af 


{C24) 


1  -  u*  =  f 


(C25) 


and  have 


q-  f  —  q-  T 

w  -wFM 


^  O 


f  I  £1  e-’^A  £  ,  ,-AtA 


(C26) 


5^+ 


As  was  done  with  the  convective  heat  transfer  bridging, we  evaluate 
fiY/Ay)  in  the  simplest  way  (see  sketch  (e)): 


1  -  u*  =  f 


(C27) 


The  integration  of  equation  (C26)  then  yields 


Tw  -  't^fM  1  At 


1  -  e 


(C28) 


The  ratio,  A/A^,  becomes  very  large  in  the  free -molecule  limit  and  very  small 
in  the  continuum  limit  in  a  manner  similar  to  A/Ag^  (as  described  above  in  the 
Convective  Heat  Transfer  section).  When  A//\.  is  large,  in  equa¬ 
tion  (C28)  approaches  as  it  should.  When  A/A^  is  small,  mast 

approach  equation  (C28)  we  get 


At  t* 


(C29) 


Since  the  characteristic  thickness.  At,  has  not  been  specified,  it  can  now  be 
defined  as  having  values  that  satisfy  equation  (C29) •  We  are  thus  forcing  our 
bridging  relationship  to  have  the  correct  continuum  asymptote. 

We  can  show  directly  that  equation  (C29)  is  a  consistent  relationship  by 
making  use  of  equations  (5-6,  6.23,  6.25,  6.29,  6.30,  6.32,  and  6.4l)  of 
reference  30.  We  can  also  show  the  consistency  of  equation  (C29)  by  using  the 
heat-transfer  relations  that  we  have.  By  combining  equations  (l6) ,  (17),  (35), 
and  (36)  we  obtain 


(C30) 


AsAe^l  1  + 


A-cm  V  P21 


Then,  using  equation  (C19b) ,  we  have 


•^sA-ca  {  ^ 


■cm 


(C31) 


This  gives  equation  (C29)  when 


55 


^cm 


n/  P21 


A3A, 


cq. 


Ea 

1  +  ^  t 


(C32) 


When  we  substitute  equation  (C29)  into  (C28) ,  we  have  the  bridging 
relation  for  the  x  gradient  of  surface  shear. 


T  *  =  T  * 

'w  Vc 


1  -  e 


-TwFm/' 


wc 


(37) 


If  functional  forms  for  fiY/fSry)  other  than  the  simple  one  in  equation  (C27) 
are  chosen  for  the  derivation  given  above ^  the  bridging  equations  will  be  more 
complicated  than  equation  (37)  but  similarly  behaved.  It  is  possible  that 
some  of  the  more  complicated  forms  may  furnish  more  accurate  bridgings  for 
some  specific  shear  situations. 

Direct  experimental  verifications  of  equation  (37)  applied  to  the 
stagnation  region  with  blowing  have  not  been  found.  However^  it  is  believed 
that  this  relationship  is  a  reasonable  and  useful  approximation  for  surface 
shear  in  the  transitional  regime.  The  expression  is  a  monotonic  function  of 
t^Pm/'^wc  correct  free -molecule  and  continuum  regime  asymptotes. 

Although  equation  (37)  ^sis  derived  for  a  stagnation  region  boundary 
layer it  is  of  possible  interest  to  compare  calculated  results  with  subsonic 
flat  plate  shear  measurements  since  these  data  are  available.  This  comparison 
is  shown  in  figure  8.  The  measured  data  are  from  figure  9(a)  of  reference  32. 

The  agreement  shown  in  figure  8  is 
thought  to  be  surprisingly  good. 

An  attempt  was  also  made  to 
compare  calculated  shear  from  equa¬ 
tion  (37)  "With  measured  shear  in  low 
speed  Couette  flow  as  reported  in 
reference  33  (not  shown) ;  the  agreement 
obtained  was  fair.  An  equation  derived 
for  Couette  flow  by  the  authors  of 
reference  33  checks  the  data  very 
closely;  this  illustrates  that  one 
should  generally  prefer  a  bridging 
equation  derived  for  a  specific  situa¬ 
tion  if  this  be  available. 


Figure  8.-  Comparison  of  calculated  skin  friction 
bridging  with  measured  skin  friction  for 
subsonic  flat  plate  (without  blowing) . 


Drag  Bridging 

For  flight  calculations^  we  use  the  trajectory  equations  (52)  to  (5^)^ 
in  which  the  quantity^  M/CpA^  appears.  The  evaluations  of  Cp  and  M/CpA  are 
given  in  equations  (57)  Q^nd  {59)?  respectively^  equation  (57)  being  the 
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■bridging  formla  for  Cj).  A  'bridging  relation  for  the  drag  coefficient  is 
clearly  necessary  for  entry  flights  since  Cj)  mst  initially  have  a  free- 
molecule  evaluation,  and,  for  many  entries,  the  final  value  of  Cj)  -will  he 
the  continuum  value. 

To  evaluate  the  drag  bridging,  we  again  make  use  of  our  first-collision 
model  (in  a  treatment  that  amounts  to  a  further  approximation  with  the  model). 
We  again  use  an  unspecified  characteristic  effective  collision  thickness,  Ap, 
and  assume  that  collisions  outside  of  Ap  have  a  negligible  effect  on  drag. 

We  assume,  then,  that  the  molecules  that  have  a  longer  free  path  than  Ap 
make  their  first  effective  collision  with  the  body  and  contribute  to  free- 
molecule  drag.  The  other  molecules  collide  with  each  other  within  the  depth, 
these  molecules  bathe  the  body  in  a  continuum  fluid  and  contribute  to  con¬ 
tinuum  drag.  A  more  rigorous  development  would  require  Ap  to  vary  with  posi¬ 
tion  on  the  body,  but  we  will  take  Ap  to  be  some  average  value  for  the  whole 
body.  Similarly,  the  value  of  the  mean  free  path  of  the  molecules  near  the 
body  will  depend  on  position  on  the  body,  but  we  will  use  a  nominal  or 
averaged  mean  free  path,  Aa,  for  the  gas  between  the  shock  wave  and  the,  body. 


According  to  our  model,  we  can  sum  the  free -molecule  and  continuum  drags 
and  obtain 

(C33a) 


Cp  =  CpFM  (  )  +  Cpc  J 


(C33b) 


The  dynamic  pressure  ratios  are  evaluated  as  density  or  mass  fraction  ratios. 

SxFM  _  ^  -Ap/Aa 

^  ^00 


^  =  1  .  (C3hb) 

Qoo  ^00 


Combining  eq[uations  (0314-)  witb  (C33b)^  we  have 

We  can  insert  the  quantity.  Eg,  from  equation  (58)  into  equation  (C35)  to 
obtain 


(C35) 


Cp  =  Cpc  (1  +  Ege-^C/^^^ 


(036) 
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We  observe  that  equations  (C35)  aJ^d-  (C36)  have  the  correct  asymptotes.  The 
variable,  is  a  reciprocal  Knudsen  number.  To  this  point,  the  deriva¬ 

tion  is  similar  to  that  found  in  reference  3^j  "bu-t  from  a  somewhat  different 
point  of  view.  In  reference  34,  Ap  is  specified  as  the  shock  standoff  dis¬ 
tance,  evaluated  empirically  in  terms  of  and  R  for  a  sphere;  this  recip¬ 

rocal  Knudsen  number  is  then  multiplied  by  a  constant  factor  to  fit 
experimental  data. 

We  can  postulate  a  functional  form  for  Ap  as 

^  ^  'body  shape) 


From  dimensional  considerations^  ve  write 


Ad  =  "body  shape) 


where  the  nose  radius,  R,  characterizes  the  body  size.  Then  we  have 


We  can  write 


^  ^  R_ 
As  As 


f 


P 


As  As  ^2 


(C37a) 


(C37b) 


This  latter  form  is  equivalent  to  the  exponent  in  equation  (7)  of  reference  34 
when  p  is  taken  as  a  constant  for  a  spherical  body  shape .  As  shown  in 

reference  34,  with  the  constant  properly  adjusted,  a  good  fit  is  obtained  with 
the  drag  data  for  a  sphere  at  high  speeds  in  air  and  in  helium. 


For  practical  calculations,  it  seems  preferable  to  avoid  evaluating  A^. 
We  can  use  the  classical  relation  between  viscosity  and  mean  free  path 
(eq.  (119)  of  ref.  27)  which  states  that  viscosity  is  proportional  to  the 
product  of  Ap  and  a  mean  molecular  speed.  We  can  express  A2P2  in  terms  of 
an  undissociated  value,  A^^p^^;  this  involves  using  fpCpts,  Tg),  the  ratio  of 
actual  viscosity  to  the  undissociated  (Sutherland)  viscosity.  (Values  of  fn 
for  air  are  tabulated  in  table  VI  of  ref.  350  We  make  use  of  the  Sutherland 
formula. 


i!-  = 


(C38) 


to  relate  to  free -stream  conditions^  and  we  obtain  the  approximation 


^2^2  J 


Ts  +  b(Tg/T^) 
Ts  +  b 


(C39) 
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Using  equation  (l06)  of  reference  27 ^  ve  can  write 


.  - 
^Poo  q 


(c4o) 


and  put  equation  (C39)  in  the  form 


+  h(Ts/Too) 


^2^2  "  n/  “^co 


Ts  +  b 


(Chi) 


Substituting  this  expression  into  equation  (C37b)  we  have 


(RD)lO~%/p 


^2  Voo  fTs  +  ^(T^/^oo) 


So3  n/  lHa  To  +  b 


(C42) 


We  let 


_ IQ-^Pgifp 

f'J^s  +  ^(^sAco) 


^00  N  ^  L  “^S  ^ 


(Clf3) 


and  we  have 


Aj) 


ickk) 


where 


fy  =  fy(Reoo,  Mco^  gas,  body  shape) 


(C45) 


For  given  flight  conditions  and  gas,  f^  will  be  a  function  of  body  shape 
(includes  angle  of  attack).  For  many  flight  cases,  body  shape  does  not  change 
mach  in  the  passage  through  the  transitional  regime;  it  can  be  expected  that 
generally  fy  -will  not  undergo  a  large  variation  in  the  transitional  passage. 
We  now  approximate  fy  as  a  constant  for  the  flight  of  a  given  body  and  use 


fy  =  15(1  +  E14) 


(C46) 


where  E14  can  be  assigned  the  value  zero  for  a  sphere  in  air,  and  the  con¬ 
stant,  15,  has  been  selected  to  match  experimental  data  as  described  below. 
We  combine  equations  (G36),  (C44),  and  (C46)  to  obtain  our  drag  bridging 
equation. 
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The  parameter^  E14,  can  he  adjusted^  in  principle^  to  account  approximately 
for  the  effects  of  all  the  variables  in  equation  (C45)  (during  passage  through 
the  transitional  regime).  Unless  data  are  at  hand^  E14  for  air  would  nor¬ 
mally  he  assigned  the  value  zero^  as  for  a  sphere.  Although  the  (presumably 
modest)  dependence  of  E14  on  body  shape  may  not  be  known^  the  values  of  Cp 
calculated  from  equation  (5?)  will  also  depend  directly  on  body  shape  through 
CpQ  and  E9.  The  constant^  15^  in  equations  (C46)  and  (5?)  'W'as  determined  by 
matching  calculated  values  of  drag  with  the  measured  drag  data  for  spheres  in 
air  of  reference  36*  The  comparison  is  shown  in  figure  9-  Corresponding  with 

calculated  values  of  Cp  for  figure  9^ 
the  values  of  RD  used  in  equa- 
tlon  (57)  were  converted  to  Reoo^ 
because  the  data  of  reference  36  are 
plotted  with  Reoo  as  abscissa.  The 
data  of  reference  36  were  obtained  in 
air  at  a  nominal  settling  chamber  tem¬ 
perature  of  300^^  K  over  a  Mach  number 
range  from  3*8  to  4.3*  The  data  indi¬ 
cate  that  any  Mach  nimiber  effect  on 
the  transition  of  the  drag  coefficient 
is  probably  small.  The  value^  E14  =  0 
was  also  used  to  obtain  a  good  fit  with 
sphere  drag  data  in  air  reported  in 
reference  34  (not  shown) .  One  set  of 
these  data  was  obtained  in  undissociated  air  at  a  nominal  settling  chamber 
temperature  of  2500^  K  over  a  Mach  nxmiber  range  from  I5.98  to  20. 90.  Another 
set  of  data  was  obtained  in  dissociated  air  over  a  range  of  hypersonic  Mach 
numbers  (11.34  to  58*7)  at  a  nominal  settling  chamber  temperature  of  9000*^  K. 

A  third  set  of  sphere  drag  data  was  obtained  in  helium  with  hypersonic  flow 
conditions.  The  helium  data  require  E14  2  for  a  good  fit. 

The  drag  coefficient^  Cp^  is  used  only  in  equation  (59)  for  the  quantity 
M/CpA.  Equation  (59)  contains  an  empirical  evaluation  of  the  variation  of 
body  mass^  M,  with  the  stagnation  point  surface  recession^  X.  Any  change  in 
the  continuum  drag  coefficient^  Cpc^  can  also  be  accounted  for  in  the 
empiricism  of  the  M/CpA  evaluation. 


Figure  9. -  Comparison  of  calculated  and  measured 
values  of  the  drag  coefficient  for  spheres  in 
the  transitional  regime. 


Example  of  Free -Molecule -Continuum  Bridging 

The  bridging  relations  developed  in  this  appendix  are  essentially 
approximations  to  be  used  for  the  conditions  for  which  they  were  derived. 

For  other  conditions^  more  rigorously  developed  equations  may  be  available^ 
as,  for  example,  for  cases  of  shear  in  Couette  flow  such  as  reported  in  refer¬ 
ence  33*  Under  conditions  for  which  no  well  developed  bridging  relations 
exist,  the  formulas  given  in  this  appendix  are  recommended  as  engineering 
approximations.  The  structure  of  these  equations  insures  that  the  equations 
have  the  correct  asymptotes,  and  this  tends  to  limit  inaccuracies. 
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DRAG,-^ 

COFM 

-HEAT 

JRANSFER, 


VAPORIZATION;^ 


Representative  curves  obtained 
BLOCKAGE  FACTOR  j  froin  thc  bridging  ec[uations  are  pre- 

:  sented  in  figure  10.  The  ordinate 

u  '  quantities^  normalized  to  their  free- 

<  -  Jsfer  molecule  values^  are  plotted  against  a 

t  ^  —  reciprocal  Knudsen  number  which  was 

ia‘'r  varied  by  varying  the  nose  radius^  R. 

°  i  -  vaporizationXs^s^  Although  the  bridging  equations  are 

i  I  general,  the  numerical  inputs  to  the 

equations,  and  therefore  the  curves 
obtained,  will  depend  on  the  individual 

01 1-  I  I  '  '  '  ^ 

■  -01  -I  1-0  10  100  1000  case  calculated.  The  calculated  curves 

BODY  RADIUS  „ _  _  _  .  .  „ 

MEAN  FREE  PATH  ’  lu  the  flgure  are  for  a  tektite  glass. 

^  114--  4>  A  clirve  for  the  quantity  q,,-.  (eq.  (l8)) 

"bridging  between  free -molecule  and  iS  not  Shownj  this  quantity _ iS  the 

continuum  flow.  product  of  the  quantities,  ^  and  q^^ 

(see  eq.  {22b)),  for  which  the  curves  are  shown.  The  dotted  lines  on  the 
right  side  of  the  figure  are  the  continuum  regime  asymptotes  for  the  various 
quantities,  while  on  the  left  side  the  free-molecule  asymptote  for  all  quanti¬ 
ties  is  unity.  It  is  noted  that  the  various  quantities  approach  their  asymp¬ 
totes  at  different  rates.  (Under  some  conditions,  one  can  consider  that  all  of 
the  quantities  plotted  are  not  even  in  the  same  regime.)  The  figure  illus¬ 
trates  probably  the  most  important  feature  of  the  bridging  relations  used. 

They  automatically  place  control  of  the  various  quantities  in  the  appropriate 
controlling  regime,  the  free-molecule,  transitional,  or  continuum. 


1.0  10 
BODY  RADIUS  ^ 
MEAN  FREE  PATH  ’ 


Figure  10.-  Representative  calculation  of 
bridging  between  free-molecule  and 
continuum  flow. 


Curves  corresponding  to  those  of  figure  10  can  be  plotted  for  other 
materials,  and  most  of  the  curves  will  be  similar,  but  somewhat  displaced. 

The  heat  transfer  curve  (q^^)  is  approximately  universal;  the  drag  curve  (Cp) 
is  influenced  considerably  by  body  shape;  the  other  curves  (t^,  v^,  ^Ir)  will 
vary  somewhat  with  relative  rates  of  vaporization  (or  reaction)  of  the 
material  being  considered. 
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APPENDIX  D 


USE  OF  COMPUTING  PROGRAM 


The  computing  program  can  be  used  to  solve  a  variety  of  problems 
involving  surface -type  ablation.  Data  are  arbitrarily  read  in  subject  to  the 
limitations  given  below  in  the  listing  of  Input  Data.  As  explained  in  the 
section  METHOD  OF  THE  NUMERICAL  PROGRAM^  the  stability  parameter,  Z,  should  be 
<  1/2  for  all  grid  points  and  at  all  times.  The  quantity,  Z,  is  printed  out 
by  the  program  so  that  its  values  may  be  observed. 


Computing  Program  Options 

The  numerical  computing  program  has  six  major  groups  of  options  as  listed 
below. 

1.  Running  conditions 

(a)  Normal  wind  tunnel,  KF  =  1. 

(b)  Rarefied  wind  tunnel,  KF  =  2. 

(includes  all  wind  tunnel  cases,  but  computing  time  is  longer  than 
with  option  (a)) 

(c)  Flight,  KF  =  3- 

Both  of  the  wind-tunnel  options,  (a)  and  (b),  allow  wind-tunnel  condi¬ 
tions  to  be  changed  once,  if  desired,  and  also  for  the  wind  tunnel  to  be  shut 
off.  At  time,  “the  free-stream  density,  enthalpy  velocity,  and 

free-stream  velocity  are  changed,  respectively,  to  Dcha  ^cha 

V^cha  (OKBAR).  At  time,  toff  (OMCO),  the  wind  tunnel  is  shut  off,  and  the 
calculations  are  continued  while  the  model  cools.  If  these  changes  are  not 
desired,  the  values  of  and  t^^^  can  be  set  larger  than  the  time  corre¬ 

sponding  to  the  final  time  line  number,  NF  (see  Time  Sketch  below). 

2.  Internal  radiation 

(a)  Transparent,  KG  =  1. 

(b)  Opaque,  KG  =  2. 

(c)  Semitransparent,  KG  =  3* 

3 .  Surface  conditions 

(a)  Evaporation  or  sublimation,  KCH  =1. 

(b)  Surface  chemical  reaction,  KCH  =  2. 
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4.  Initial  conditions  for  flight  case  (KF  =  3) ;  (see  appendix  B) . 

(a)  Normal  entry  into  an  atmosphere  using  a  computed  exponential  tempera¬ 
ture  profile  in  the  hody.  With  an  assumed  Tq)^  the  initial 

values,  Di  and  Ai  are  computed  hy  the  program,  KDEL  =  1. 

(h)  Arbitrary  initial  values  of  atmospheric  density,  D^,  and  thermal 
thickness  of  exponential  temperature  profile  in  the  body,  A^, 

KDEL  =  2. 

5.  Back  face  boundary  conditions 

(a)  Back  face  aerodynamically  exposed,  KBAK  =  1. 

(b)  Backing  material  forming  a  heat  sink,  KBAK  =  2.  (For  a  heat  sink  of 

zero  heat  capacity,  or  an  adiabatic  back  boundary,  KBAK  =  1  should 
be  used. ) 

6.  Planet  and  atmosphere  for  flight  case  (KF  =  3) 

(a)  Earth  entry  with  the  ARDC  atmosphere  (approximated  exponentially  with 

3  programmed  values  of  scale  height) .  Earth  radius  programmed  at 
6440  km  (Rp  in  eq..  (53))j  KC5  =  1- 

(b)  Arbitrary  planet  with  exponential  atmosphere  having  arbitrary  scale 

height  (initial,  two  intermediate,  and  final  values),  KC5  =  2. 


Nomenclature  of  Computing  Program 

The  nomenclature  used  in  the  computing  program  is  in  symbolic  FORTRAN 
language.  Separate  listings  of  input  and  output  data  are  shown  below. 


Input  Data 

Input  data  are  listed  below  in  their  order  of  card  punching.  Actual 
card  formats  are  shown  in  the  Input  Card  Format  Sketch.  All  input  data  are 
printed  out  by  the  program  in  an  initial  readout  (see  Sample  Case,  below). 

A  quantity  listed  as  an  option  is  defined  in  the  section.  Computing  Program 
Options,  above. 

Following  the  definition  of  a  quantity,  the  value  of  an  option  selection 
may  be  shown.  The  particular  quantity  is  not  needed  (and  not  used)  for  other 
values  of  the  option  selection.  Unused  quantities  are  normally  assigned  the 
value  zero;  in  any  case,  the  input  card  formats  must  be  maintained. 

The  maximum  number  of  grid  points  to  be  used  in  the  finite  difference 
spacing  is  98^  this  is  the  maximum  value  for  the  quantity,  MF  (see  Spacing 
Sketch) . 
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CARD 


a  (6  ao  24  28  5^  36|  40|  44|  481  52|  56|  60|  64|  681  72| 


KF  KG  Nl  N2  MF  kg  K3  K4  LN  M2  M3  MF  Jl  I  J2  |KC5|KCTjKM  l|KMg 
M3  KMI  ltCII2  KCII3  KCNF  (Nl  M!  KCKI  KCK2  KCIIF  ICH  KOEl  XBAK 


ALL  NUMBERS  (N  14  FORMAT  (MUST  BE  RIGHT  JUSTIFIED) 


1 

10 

(9 

28 

37 

46 

55 

64  72 

Al 

A2 

A3 

A4 

81 

B2 

B3 

B4 

B5 

86 

87 

B8 

89 

BIO 

BM 

BIZ 

BI3 

BI4 

BIS 

BI6 

Cl 

C2 

C3 

C4 

C5 

C6 

C7 

ce 

El 

E2 

E3 

E4 

E5 

E6 

E7 

E8 

E9 

ElO 

Ell 

EI2 

EI3 

EI4 

EI5 

EI6 

SIGMA 

RO 

OMCO 

OK  BAR 

CHM 

VCINF) 

GAMAI 

TWI 

DELTW 

ALLOW 

HS 

DCI 

DTT2 

DTT  3 

OYI 

RH02I 

OMO 

EI7 

EIB 

DELTJ 

EI9 

R8 

E20 

E2I 

E22 

E23 

E24 

E25 

E26 

E27 

EZ8 

£29 

E30 

E3I 

E32 

E33 

E34 

E35 

E36 

E37 

E3S 

E39(0PE)I 

£40  (IPEIi: 

ALL  NUME 

EPS  IN 

10 

9.3  FORM 

19 

AT  (DEC 

28  36 

MAt 

40 

POINT  NEEDED) 

72 

ALPHA 

ALPHA2 

EMAX 

RN 

IT 

{E9.3) 

(E9.3) 

(E9.3) 

(E9.3) 

(14) 

FORMATS  ARE  INDICATED  BY  PARENTHESES.  E  FORMATS 
REQUIRE  DECIMAL  POINT,  I  FORMAT  MUST  BE  RIGHT  JUSTIFIED. 

Input  Card  Format  Sketch 

Ml,  M2.  M3  ARE  BREAK  POINTS  FOR  CHANGING  At] 

FINE  SPACED 
e.g.,LN=4 


n+lll 


GRID  POINT  NO.=M 


:al 

e.g.,Mi  =  i7 


(A7]J 


(A^4} 


>0^-4 

^  (CALCULATED) 


M3  MF 

(MAX  =  98) 


INITIALLY  AYj  =  At7| 

At]^,  At]^,  Arj^  ARE  INCREMENTS  USED  IN  FINITE  DIFFERENCE  EQUATION 

AT],  INCREMENT  USED  FOR  INTERPOLATIONS 

INITIAL  VALUE  OF  LENGTH/A 7],  <  1665  REQUIRED  FOR 
TRANSPARENT  CASE  (KG  =  I) 

Spacing  Sketch 


Nl,  N2  ARE  BREAK  POINTS  FOR  CHANGING  At 


- At  I  =  DDT  I  ■ 


At9=DDT2 


TIME  LINE  NQ=N  I 
t  =  ti 


At,,  Atg,  Atj,  ARE  INCREMENTS  USED  IN  FINITE  DIFFERENCE  EQUATION 

Time  Sketch 


ONE  TIME  LINE 

KMI,  KM2,  KM3  ARE  BREAK  POINTS  FOR  PRINTING  INCREMENTS 


(KCMI) 

(KCM2) 

(KCM3) 

(KCMF) 

GRID  POINT  NO.  I  KMI  KM2  KM3 

KCMI,  KCM2,  KCM3,  KCMF  ARE  GRID-POINT  PRINTING  INCREMENTS 


BETWEEN  TIME  LINES 

KNI,  KN2  ARE  BREAK  POINTS  FOR  PRINTING  INCREMENTS 


(KCNI) 

(KCN2) 

(KCNF) 

TIME  LINE  NO.  I  KNI  KN2 

KCNI,  KCN2.  KCNF  ARE  TIME  LINE  PRINTING  INCREMENTS 

Printing  Sketch 


Card  A 

(All  numbers  are  integers  in  Ik  FORMAT) 

KF  Running  condition  option 

KG  Internal  radiation  option 

Nl  Time  line  number  at  which  the  finite  time  Increment  At  (DTT)  changes 
from  Ati  (dDTI)  to  Ata  (DDT2)  (see  Time  Sketch). 

N2  Time  line  number  at  which  the  finite  time  increment  At  (DTT)  changes 
from  Ata  (DDT2)  to  Ats  (DDT3)  (see  Time  Sketch). 

NF  Final  time  line  number  (see  Time  Sketch). 

K2  Defined  by  Ari^  =  (K2)  Ati^  (or  Ay^  =  (K2)  Ay^) .  Can  be  1  for  opaque 

and  semitransparent  cases  (KG  =  2,  3).>  kiut  must  be  at  least  2  and  even 
for  transparent  case  (KG  =  l)  (see  Spacing  Sketch). 
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K3 

Ki^ 

LN 

M2 

M3 

MF 

J1 

J2 

KC5 

KCT 

KMl 

KM2 


Defined  Dy  Arjg  =  (K3)  Ati^  (or  =  (K3)  Ay^) .  Can  be  1  (see 

Spacing  Sketch) . 

Defined  by  At^^  =  (K4)  Ar\^  (or  Ay^  =  (Ki^-)  Aj^).  Can  be  1  (see 
Spacing  Sketch) . 

Increments  of  Atj^  spacing  over  which  At)^  spacing  exists. 

Must  be  >  k  (see  Spacing  Sketch). 

Grid  point  at  which  space  increment  changes  from  Ati^  to  Arig 
( see  Spacing  Sketch) . 

Grid  point  at  which  space  increment  changes  from  Aiig  to  At^^ 

(see  Spacing  Sketch) . 

Grid  point  at  back  face.  Maximum  value  =  9d  (see  Spacing  Sketch). 

Order  of  Interpolation  for  T  (TE)  in  the  fine  spaced  (Ari^) 
region. 

Order  of  interpolation  for  A|_i  (YDEL)  . 

Planet  and  atmosphere  option  for  flight  case  (KF  =  3) • 

Maximum  number  of  Iterations  to  determine  front  face  temperature^ 
T-fff  within  allowable  error  selected  (ALLOW) .  If  exceeded^ 
calculation  will  stop . 

Grid  point  at  which  printing  interval  on  one  time  line  changes 
from  KCMl  to  KCM2  (see  Printing  Sketch). 

Grid  point  at  which  printing  interval  on  one  time  line  changes 
from  KCM2  to  KCM3  (see  Printing  Sketch). 


J 

KM3 


Card  B 

(All  numbers  are  integers  in  1 4  FORMAT) 

Grid  point  at  which  printing  interval  on  one  time  line  changes 
from  KCM3  to  KCMF  ( see  Printing  Sketch) . 


KCM1^KCM2^  1  PT-intinc  intervals  of  grid  points  (see  Printing  Sketch). 
KCM3,KCMF  J  ^ 

ipiuie  line  number  at  which  time  line  printing  Interval  changes 
from  KCWl  to  KCW2  (see  Printing  Sketch) . 

10^2  Time  line  number  at  which  time  line  printing  interval  changes 

from  KCW2  to  KCNF  (see  Printing  Sketch). 
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KCNl 

KCNF 

KCH 

KDEL 

KBAK 

E9.3 

A1 

A2 

A3 

A4 

B5 

b6 

B7 

b8 

BI3 

B14 

B15 

B16 

Cl 

C2 

C3 

C4 
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j-  Time  line  printing  intervals  (see  Printing  Sketch). 
Surface  conciition  option. 

Initial  condition  option  for  flight  case  (KP  =  3). 
Back  face  Boundary  condition  option. 


For  cards  1-11^  there  are  normally  8  numbers  per  card^  each  number  in 
FORMAT. 


Card  1 


Ai^  equation  (27) 

B1 

P 

As,  equations  (lO)^(ll) 

B2 

Bs,  equation  (84) 

As,  equations  (32), (3 5) 

B3 

Bs,  equation  (84) 

A4,  equation  (15) 

B4 

B4,  equation  (85) 

Card  2 

equation  (85) 

B9 

B^y  equation  (87) 

Be,  equation  (86) 

BIO 

equation  (29) 

equation  (86) 

Bll 

Bll,  equations  (28)^  (30) 

T 

^0 

B12 

Card  3 

for  opaque  case  (KG  =  2) ;  ^xaax  semitransparent  case  (KG  =  3) 
Bi4^  equation  (85) 

for  opaque  and  semitransparent  cases  (KG  =2^  3) 

^16^  equations  (46)^  (47)^  (51);  Bi©  =  1  for  transparent  case  (KG  =  ij^ 
Bi6  =  0  for  opaque  and  semitransparent  cases  (KG  =  2^  3) 

Ci^  equation  (6o) 

equation  (6o) 

C3,  equations  (56)^  (59) 
equations  (56)^  (59) 


Card  4 


C5 

Shij,  for  flight  case  (KF  =  3)^  and  (KC5 

=  2). 

C6 

Oey  equation  (l5) 

C7 

L/Dr^  equation  (53)  for  flight  case  (KF 
cases  (KF  =  1^  2) . 

-  3) '>  Dcha  wind  tunnel 

C8 

g  /lO^^  equation  (53h)  for  flight  case 
^tunnel  cases  (KF  =  1^  2) . 

(KF  =  3);  for  vind 

El 

El,  equation  (86) 

E2 

equation  (87) 

E3 

E3,  equation  (88) 

EU 

E4,  equation  (44) 

Card  5 

E5 

E5,  equation  (44)  E9 

E9,  equation  (58) 

e6 

Eey  equation  (44)  ElO 

Eio^  equation  (88) 

E7 

Ey,  equation  (26)  Ell 

Ell,  equation  (61) 

E8 

Es^  equations  (34)^  (35)  E12 

Ei2^  equation  (6I) 

Card  6 

E13 

Ei3,  equation  (50) 

E14 

Ei4^  equations  (57)^  (59) 

E15 

Ei5,  equation  ( 10a) 

ei6 

Ei6^  equation  (l2) 

SIGMA 

a  =  1.369x10“^^ 

RO 

Rj_,  equation  (60) 

OMCO 

(M/CDQA)j_,  equation  (59)  for  flight  case  (KF  =  3);  t^^f  for  wind 
tunnel  cases  (KF  =  1,  2) . 

OKBAR 

K,  equation  (l2)  for  flight  case  (KF  = 

3)^  ’'^oocha  tunnel 

cases  (KF  =  2)  . 
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Card  7 


CHIl  Xx,  equation  (l2) 

VCIIFI  for  flight  case  (KF  =  3)j  for  vind  tunnel  cases  (KF  =  1,  2), 

Y^  =  VCINFI  until  t  >  t^ha  or  t  >  toff. 

GAMAI  7j_  for  flight  case  (KF  =  3)j  ^cha  wind  tunnel  cases  (KF  =  1,  2). 

TWI  >  Tq  for  entry  flight  case  (KF  =3^  KDEL  =  l) . 

DELTW  Initial  guess  for  incrementing  T^j_  in  iterating  for  Tw  in  the 
second  time  line  (can  be  zero). 

ALLOW 

HS  hg  for  wind  tunnel  cases  (KE  =  1,  2) ;  constant  at  read -in  value 

until  t  >  or  t  >  toff. 

DCI  Dj_  for  flight  case  (KF  =  3)j  read  in  only  when  KDEL  =  2;  with 

KDEL  =  ly  is  calculated  by  the  program.  For  wind  tunnel  cases 
(KF  =  1^  2)^  read  in;  D  =  DCI  until  t  >  or  t  > 

Card  8 

DDTl  Ati  (see  Time  Sketch)  . 

DDT2  Ats  (see  Time  Sketch). 

DI)T3  Ats  ( see  Time  Sketch)  . 

DYl  At]^  =  initial  value  of  Ay^  (see  Spacing  Sketch);  initial  value  of 

length/Ar]^  <  I665  is  required  for  the  transparent  case  (KG  =  l) . 

RH021  P21;  read  in  for  wind  tunnel  cases  (KF  =  1^  2) ;  for  flight  case 

(KF  =  3)^  P21  calculated  continuously  (eq.  (61)). 

OMO  %/F^  equation  (48);  used  only  when  back  face  is  aero  dynamic  ally 

exposed  (KBAK  =  l) . 

EI7  equation  (60) 

EI8  Eis^  equations  (56)^(59)^  used  for  flight  case  (KF  =  3). 

Card  9 

DELTJ  Initial  guess  for  incrementing  T.^  in  iteration  for  time  lines  3? 

4^  5j  should  not  be  zero. 

EI9  equation  (49)  j  used  when  back  face  is  aerodynamically  exposed 

(KBAK  =  1) . 
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KB  Rb^  equations  (ijS),(50)]  used  when  back  face  is  aerodynamically  exposed 

(KBAK  =  1) . 

E20  Eso^  equation  (87) 

E21  equation  (ij),  (KF  =  2,  3);  if  unknown  use  A(,q  =  1. 

E22  Acm^  equation  (36),  (KP  =2,  3);  if  unknown  use  Acm  =  1- 

E23  A^y)  equations  (39c,  39d),  (KCH  =1);  if  unknown,  evaluate  using 
Aqy  =  1  in  equation  (39c) . 

E24  B,  equation  (39e),  (KCH  =  2). 

Card  10 

E25  Eij,  equation  (39®) ,  (KCH  =  2)  . 

E26  c",  equation  (51)^  (KBAK  =2);  (if  c’=0,  use  KBAK  =  1  and  OMO  =  O) . 

E27  Rp,  equation  (53) ^  flight  case  (KF  =3)^  and  (KC5  =  2). 

E28  flight  case  (KF  =  3)^  and  (KC5  =  2). 

E29  flight  case  (KF  =  3)^  and  (KC5  =2). 

E30  Sbg,  flight  case  (KF  =  3)^  and  (KC5  =  2). 

E3I  ,  flight  case  (KF  =  3)^  and  (KC5  =  2). 

“*■0 

E32  in  equation  (BIO)^  flight  case  (KF  =  3);  used  only  "v^hen  KDEL  =  1 

and  KC5  =  2;  has  arbitrary  value^  hut  may  equal  (C5)  • 

Card  11 

E33  E33,  equation  (45)^  semitransparent  case  (KG  =3)* 

E34  in  equation  (Bl)  ;  read  in  only  for  flight  case  (KF  =  3)^  with 

arbitrary  initial  conditions  (KBEL  =  2) .  Otherwise  the  program 
computes  from  equation  (B2b) . 

E35  E35,  equation  (28) 

E36  E36^  transparent  case  (KG  =  l) ;  equation  (97) i  E36  =  1*0  has  given 

good  energy  balances. 

E37  ^37^  transparent  case  (KG  =  l) ;  equation  (98);  E37  =  1.0  has  given 

good  energy  balances. 

E38  E38^  flight  case  (KF  =  3);  equation  (55)- 

E39^E1}-0  (Open) 
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Card  12 


This  card  is  only  required  for  the  transparent  case  (KG  =  l) .  The  format 
for  this  card  is  shown  in  the  Input  Card  Format  Sketch.  The  FORTRAN  quanti¬ 
ties  ALPHA^  ALPHAS,  EMAX,  RM  are  in  FORMAT  E9.3;  IT  is  in  FORMAT  Ik. 

ALPHA  a,  equation  (4l) 

ALPHAS  equation  (4l) 

EMAX  Smax^  equation  (4S) 

RN  n,  equations  (41,  4S) 

IT  Spacing  of  calculation  and  printing  of  F  and  g  =  hF/hy  in  the 

fine  spaced  region  ( see  Spacing  Sketch) .  IT  must  he  an  integer 
factor  of  KS  hut  must  he  less  than  KS;  for  example,  if  KS  =  6, 

IT  can  he  1,  S,  or  3-  IT  =  3  means  that  F  and  g  are  printed 
out  for  every  third  grid  point  in  the  fine  spaced  region;  in  the 
remaining  regions  they  are  printed  for  every  grid  point.  The 
value  of  IT  affects  somewhat  the  accuracy  of  the  finite  differ¬ 
ence  and  energy  balance  calculations.  IT  =  1  gives  the  greatest 
accuracy.  A  larger  value  of  IT  decreases  accuracy  while  reducing 
computing  time. 


Output  Data 

All  input  data  are  printed  out  in  an  initial  output  format  (see  Sample 

Case  he  low) .  In  addition,  the  following  quantities  (some  calculated)  are 

printed  in  an  initial  output  (see  Sample  Case  helow) . 

QOII  q^i 

Q02I  qj^^ 

TvRI 

111  t.;  equation  (b4)  for  wind  tunnel  cases  (KF  =  1,  2);  equals  zero  , 

^for  flight  case  (KF  =3)* 

DELI  A;j^;  equation  (B2h)  for  wind  tunnel  cases  (KF  =  1,  2)  and  flight 

case  (KF  =  3)  with  KDEL  =  1;  for  flight  case  with  KDEL  =  2,  Aq 
equals  E34. 

DCI  Dq;  read  in  quantity  for  wind  tunnel  cases  (KF  =  1,  2)  and  flight 

case  (KF  =  3)  with  KDEL  =  2;  from  equation  (BIO)  for  flight  case 
with  KDEL  =  1. 

RHO  =  D/1226  (for  this  printing  D  =  Dq). 

Ml  Ml  =  1  +  (K2) (LN) ;  see  Spacing  Sketch. 
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The  spacing  of  printing  intervals  is  arbitrarily  selected  (see  subsection 
Input  Data  and  the  Printing  Sketch) .  For  each  time  line  printed  out,  the 
following  quantities  appear.  (See  section,  Sample  Case,  below  for  format.) 


N 

time  line  number 

CMC 

M/Cj)A 

TT 

t 

XIP 

%r  = 

R 

R 

cm 

X 

SRGM 

sin  7 

DTT 

At 

RHO 

Poo 

PDP 

P" 

VC 

V 

P 

P 

VCINF 

V 

00 

PHI 

9 

PSI 

V 

A2D 

■^D 

PSIBAR 

FEAR 

F 

QOC 

*^0 

FIHT 

^t 

QOFM 

ADGABS 

apg 

QO 

lo 

AGABS 

ag 

QOO 

loo 

EPSRER 

^BF 

PSQOO 

EPSLN 

^FF 
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PT2 

Pt2 
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iR 
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'^stor 
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TAUFMP 
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DTDYW 

PSQO 

FR¥ 

DVAP 

DSTOR 


^con 

^rad 

^vap 

%tOT’  t 

vhen  DSTOR  =  (stored 
energy  at  tj_)/At;  stored 
energy  at  ti  =  absolute 
value  of  second  term  in 
eq.uation  (82d) . 


^vcon 

%Lcon 

^res 

KC^  number  of  iterations 
req.uired  to  obtain  T^.. 

Tv 

1  "Vcl 

I  \fm1 


VARG 

YTNT 

DEES 

KC 

T¥ 

WC 

VWFM 

W 


Also  for  each  time  line,  the  following  arrayed  quantities  print  out  in 
columns  for  the  selected  print  spacing  of  grid  points  (see  Printing  Sketch 
and  Sample  Case  helow  for  format). 


M  M  (grid  point  number) 

Y  y 

TE  T 


UDB  u 


TBMi  V 

XRAT  Xj. 

G  g 

FR  F 

Z  Z,  must  be  <  1/2  for  stability. 


Sample  Case 

To  illustrate  the  use  of  the  computing  program,  a  sample  case  has  been 
selected  which  describes  the  entry  of  a  transparent  tektite  into  the  Earth’s 
atmosphere.  Figure  11  shows  the  input  data  for  the  sample  case)  figure  12 
shows  the  printing  out  of  the  input  data  and  the  initial  calculated  values; 
figure  13  shows  the  printing  out  for  a  typical  time  line. 
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''The  aeronautical  and  space  activities  of  the  United  States  shall  be 
conducted  so  as  to  contribute  ...  to  the  expansion  of  human  knowl¬ 
edge  of  phenomena  in  the  atmosphere  and  space.  The  Administration 
shall  provide  for  the  widest  practicable  and  appropriate  dissemination 
of  information  concerning  its  activities  a72d  the  results  thereof.^^ 

— National  Aeronautics  and  Space  Act  of  1958 
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